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Abstract The first stars to form in the Universe - the so-called Population III stars 
- bring an end to the cosmological Dark Ages, and exert an important influence 
on the formation of subsequent generations of stars and on the assembly of the first 
galaxies. Developing an understanding of how and when the first Population III stars 
formed and what their properties were is an important goal of modem astrophysical 
research. In this review, I discuss our current understanding of the physical pro- 
cesses involved in the formation of Population III stars. I show how we can identify 
the mass scale of the first dark matter halos to host Population III star formation, 
and discuss how gas undergoes gravitational collapse within these halos, eventu- 
ally reaching protostellar densities. I highlight some of the most important physical 
processes occurring during this collapse, and indicate the areas where our current 
understanding remains incomplete. Finally, I discuss in some detail the behaviour 
of the gas after the formation of the first Population III protostar. I discuss both the 
conventional picture, where the gas does not undergo further fragmentation and the 
final stellar mass is set by the interplay between protostellar accretion and protostel- 
lar feedback, and also the recently advanced picture in which the gas does fragment 
and where dynamical interactions between fragments have an important influence 
on the final distribution of stellar masses. 



Simon Glover 

Universitat Heidelberg, Zentrum fur Astronomie, Institut fur Theoretische Astrophysik, Albert- 
Ueberle-Str. 2, 69120 Heidelberg, e-mail: glover@uni-heidelberg.de 



1 



2 



Simon Glover 



1 Formation of the first star-forming minihalo 

1.1 The Jeans mass and the filter mass 

In the currently dominant ACDM paradigm, gravitationally bound objects form in 
a hierarchical fashion, with the smallest, least massive objects forming first, and 
larger objects forming later through a mixture of mergers and accretion. The mass 
scale of the least massive objects to form from dark matter is set by free-streaming 
of the dark matter particles, and so depends on the nature of these particles. How- 
ever, in most models, this minimum mass is many orders of magnitude smaller than 



the mass of even the smallest dwarf galaxies (Green, Hofmann & Schwarz 2005). 
More relevant for the formation of the first stars and galaxies is the mass scale of 
the structures (frequently referred to as dark matter 'minihalos') within which the 
baryonic component of matter, the gas, can first cool and collapse. 

A lower limit on this mass scale comes from the theory of the growth of small 



density perturbations in an expanding universe (see e.g. Barkana & Loeb 2001 1. 



From the analysis of the linearized equations of motion, one can identify a critical 
length scale, termed the Jeans length, that marks the boundary between gravita- 
tionally stable and gravitationally unstable regimes. The Jeans length is given (in 
physical units) by 



V GPo 

where c s is the sound speed in the unperturbed intergalactic medium and po is the 
cosmological background density. Perturbations on scales A > Aj are able to grow 
under the influence of their own self-gravity, while those with A < Aj are prevented 
from growing by thermal pressure. We can associate a mass scale with Aj by simply 



taking the mass within a sphere of radius Aj/2 (Barkana & Loeb 2001 1: 



An /Aj° 



This mass, termed the Jeans mass, describes the minimum mass that a perturbation 
must have in order to be gravitationally unstable. 

In the simplest version of this analysis, the value used for the sound speed in 
the equations for the Jeans length and Jeans mass is the instantaneous value; i.e. 
to determine Aj and Mj at a redshift z, we use the value of c s at that redshift. In 
this approximation, the Jeans mass is given in the high redshift limit (where the 
gas temperature is strongly coupled to the cosmic microwave background [CMB] 



temperature by Compton scattering) by the expression (Barkana & Loeb 2001 1 

'„2\-V2 



Mj = 1.35 x 10 s [ --^j ) M , (3) 
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where Q m is the dimensionless cosmological matter density parameter, and h is the 
value of the Hubble constant in units of lOOkms -1 Mpc -1 . In the low redshift limit 
(where the coupling between radiation and matter is weak and the gas temperature 
evolves adiabatically), the Jeans mass is given instead by 

M, =5.18 X Itf. ( ^Y" 1 (^\ ( !±i V' 2 „ e , (4) 

J V o.i5 J Vo.o26y v io y ° 

where £2\, is the dimensionless cosmological baryon density parameter. The evolu- 
tion of Mj with redshift is also illustrated in Figure [T] 



CD 




Redshift 



Fig. 1 Evolution with redshift of the Jeans mass (solid line), the filter mass computed in the limit 
where the relative streaming velocity between gas and dark matter is zero (dashed line) and the 
filter mass computed assuming a streaming velocity v = <T v b c (dash-dotted line). Also plotted is 
the critical minihalo mass, Merit* required for efficient tb cooling (dotted line). The estimate of 
the filter mass in the no streaming limit comes from Naoz & Barkana ( 2007), who account for a 
number of effects not treated in the original Gnedin & Hui ( 1998 1 formulation, while the estimate 
of Mp in the streaming case comes from Tseliakhovich, Barkana & Hirata (2011). The value of 
M cr i t was computed using Equation[38] 



A more careful treatment of the growth of linear density perturbations accounts 
for the fact that the sound speed, the Jeans length and potentially also the Jeans mass 
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may all change significantly during the time it takes for a perturbation to grow into 
the non-linear regime. |Gnedin & Hui ( 1998 ) showed that in this case, the appropriate 
mass scale separating the gravitationally stable and gravitationally unstable regimes 
is a form of time-averaged Jeans mass that they denote as the "filter mass", M F . This 
is given in physical units by 



4-71 / A F 

M F = T Po y 



(5) 



where the filter wavelength A F is given in the high redshift limit by ( Gnedin 2000 ) 

3 



A F 2 = 



1 



1- 



1 



1 



1/2' 



dz'. 



(6) 



It is possible to improve further on this analysis by accounting for spatial vari- 
ations in the sound speed (Barkana & Loeb 2005] Naoz & Barkana 2005 ) and by 
properly accounting for the separate rates of growth of the dark matter and baryonic 
perturbations in the high redshift limit in which the gas is mechanically coupled 
to the CMB by Compton scattering ( |Naoz & Barkana| |2007[ ). The net result is to 
somewhat lower the filter mass in comparison with the predictions of Equations [5} 
[6] Comparing the resulting filter mass with the Jeans mass (Figure 1), we see that the 
filter mass can be a factor of a few smaller than the Jeans mass at high redshift, but 
that for redshifts below z ~ 50, the filter mass is the larger of the two mass scales. 

Another complication was recently pointed out by Tseliakhovich & Hirata ( 2010| l. 
They show that prior to the recombination epoch, the strong coupling between gas 
and radiation leads to the gas developing a non-zero velocity relative to the dark 
matter. While the gas and radiation are coupled, the sound-speed in the gas is ap- 
proximately c/-\/3. where c is the speed of light, and the relative velocity between 
gas and dark matter is highly subsonic. Once the gas and radiation decouple, how- 
ever, the sound-speed of the gas decreases enormously, becoming ~ 6kms~ 1 at 



the end of the recombination epoch. Tseliakhovich & Hirata (2010) show that at the 
same time, the RMS velocity of the gas relative to the dark matter is about 30kms~ 1 , 
implying that the gas is moving supersonically with respect to the dark matter. The 
coherence length of the supersonic flow is of the order of the Silk damping scale 
( |Silk||1968| ), i.e. several comoving Mpc, and so on the much smaller scales corre- 
sponding to the formation of the first star-forming minihalos, the gas can be treated 



as being in uniform motion with respect to the dark matter. Tseliakhovich & Hirata 



( 2010| ) also show that the relative velocity between gas and dark matter acts to sup- 
press the growth of small-scale structure in both components, and that because this 
effect is formally a second-order term in cosmological perturbation theory, it was 
not included in previous studies based on linear perturbation theory. 

In a follow-up study, Tseliakhovich, Barkana & Hirata ( 20 1 1 ) > improve on the 
Tseliakhovich & Hirata ( 2010| ) analysis by accounting for spatial variations in the 



sound speed, and study the effect that the relative velocity between the gas and the 
dark matter has on the size of the filter mass. The magnitude of the relative velocity 
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v is randomly distributed with a Gaussian probability distribution function (PDF) 



with total variance (J^ hc , 



Pvbc(v) 



2 ^vbc 



3/2 



47Tv exp 



3v 2 
2cr v 2 bc 



(7) 



Tseliakhovich, Barkana & Hirata (2011 1 show that for a relative velocity v = a v b c 



(i.e. a one sigma perturbation), the effect of the relative velocity between gas and 
dark matter is to increase Mp by roughly an order of magnitude, as illustrated in 
Figure [T| Higher sigma pertur bations lead to even greater increases in Mp, but |Tseli-] 
|akhovich~ B arkana & Hirata (201 l| l show that the global average case (obtained by 
computing Mp for a range of different v and then integrating over the PDF given in 
Equation|7| is very similar to the one sigma case. Numerical studies of the effects of 
these streaming velocities (see e.g. Stacy, Bromm & Loeb 2011 Greifetal. 201 lb] > 
have generally confirmed this result, although these studies still disagree somewhat 
regarding the influence of the streaming velocities on minihalos with masses greater 
than the revised Mp. 

Nevertheless, even the most careful version of this analysis only tells us the mass 
scale of the first gravitationally bound structures to have a significant gas content, 
which is merely a lower limit on the mass scale of the first star-forming minihalos. 
The reason for this is that for stars to form within a minihalo, it is not enough that 
the gas be gravitationally bound; it must also be able to cool efficiently. In order 
for the gas within a minihalo to dissipate a large fraction of its gravitational binding 
energy - a necessary condition if pressure forces are not to halt the gravitational 
collapse of the gas (Hoyle 1953| Rees 1976 Rees & Ostriker 1977| > - it must be 
able to radiate this energy away. The timescale over which this occurs is known as 
the cooling time, and is defined as 



'cool 



1 n m kT 

y—l A ' 



(8) 



where n tot is the total number density of particles, y is the adiabatic index, k is 
Boltzmann's constant, T is the gas temperature and A is the radiative cooling rate 
per unit volume. If the cooling time of the gas is longer than the Hubble time, then it 
is very unlikely that the minihalo will survive as an isolated object for long enough 
to form stars. Instead, it is far more likely that it will undergo a major merger with 
another dark matter halo of comparable or larger mass before any of its gas has 
cooled significantly, since major mergers occur, on average, approximately once per 
Hubble time ( |Lacey & Cole| |1993| l. Therefore, to determine the minimum mass 
of a star-forming minihalo, we must first understand how cooling occurs within 
primordial gas, a topic that we explore in the next section. 
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1.2 Cooling and chemistry in primordial gas 



At high temperatures (T ~ 10 4 K and above), primordial gas can cool efficiently 
through the collisional excitation of excited electronic states of atomic hydrogen, 
atomic helium, and singly-ionized helium. However, it is relatively easy to show that 
most of the gas within a minihalo with M ~ Mp will have a temperature significantly 
below 10 4 K. If we assume that the gas within the minihalo relaxes into a state of 
virial equilibrium, such that the total potential energy W and total kinetic energy K 
are related by W = —2K, then we can use this fact to define a virial temperature for 
the minihalo ( |Barkana & Loeb||200l) 



T ■ — 



2k 



(9) 



where jj. is the mean molecular weight, m p is the proton mass, and v c is the circular 
velocity of the minihalo. This can be rewritten in terms of the redshift z and the mass 
M of the minihalo as 



r vir = i. 



x 10 4 



0.6 



M 



lO^-'M, 



2/3 



£2 m A c 



An(z) 187T 2 



1/3 



1 



10 



K, (10) 



where Q m (z) is the dimensionless cosmological density pa rameter evaluated at red - 
shift z and A c = 18k: 2 + 82d - 39d 2 , with d = Q m {z) - 1 ( |Bryan & Norman||l998| . 
In the standard A CDM cosmology, Q m (z) 



1/3 

brackets reduces to QJ . If we rearrange Equation 
of a cloud that has a virial temperature r v ; r = 10 4 K and that can therefore cool via 
atomic excitation, we find that 



1 at z > 6, and hence the term in square 
and solve for the mass M atom 



10 



Matom = 5 X 10' h 
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(JL) 

\0.6J 



-3/2 



-1/2 



1 



10 



-3/2 
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significantly larger than our estimates for Mj and Mp above. Minihalos with masses 
close to Mj or Mp will therefore have virial temperatures much less than 10 4 K, 
placing them in the regime where molecular coolants dominate. 

In primordial gas, by far the most abundant and hence most important molecule is 
molecular hydrogen, H2. The chemistry of H2 in primordial gas has been reviewed 
in a number of different studies (see e.g. Abel et al. |1997 Galli & Palla 1998 



Stancil, Lepp & Dalgarnol |1998| |Glover & Abelj |2008| >, and so we only briefly 



discuss it here. Direct formation of H2 by the radiative association of two hydrogen 

and so at low densities, most 



atoms is highly forbidden (Gould & Salpeter] |1963[ ) 
H2 forms via the reaction chain (McDowell 1961[ |Peebles & Dicke||1968] l 



H 
H 



e ->■ H" 
f H — ► H2 - 



(12) 
(13) 



with a minor fraction forming via the reaction chain (Saslaw & Zipoy 1967| > 
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H 2 + +r, 



HJ +H -> H 2 +H H 



(14) 
(15) 



In warm gas, H 2 can be destroyed by collisional dissociation (see e.g. Martin, Keogh 



& Mandy 1998 1 



or by charge transfer with H + (see e.g. 

H 2 + H+ 



H 2 +H -> H + H + H, 
H 2 +H 2 ^H + H + H 2 , 

Savin etal.| [2004] 



H+ 



-H, 



(16) 
(17) 



(18) 



but at low temperatures there are no collisional processes that can efficiently remove 
it from the gas. 

When the fractional ionization of the gas is low, the rate at which H 2 forms is 
limited primarily by the rate at which H ions form via reaction 12 as any ions that 



ions form via reaction 
form are rapidly converted to H 2 by associative detachment with atomic hydrogen 
(reaction [T3]>. If the fractional ionization is large, on the other hand, then many of 
the H ions formed by reaction 12 do not survive for long enough to form H 2 , but 
instead are destroyed by mutual neutralization with H + ions: 



H++H" ->-H + H. 



(19) 



The ratio of the rates of reactions 13 and 19 is given by A[T3lJH/^T9f 1 H+' where «h 
is the number density of atomic hydrogen, n H + is the number density of protons, 



and IfH] and Arrjgare the rate coefficients for reactions 13 and 19 respectively. Mu- 
tual neutralization therefore becomes significant whenever «h+/«h > ^T3^^T9\ Al- 
though the value of Arj3]/Arx9] is temperature dependent, the temperature dependence 
is weak if one uses the best available determinations of the rate coefficient s ( Kreckel| 



et al.|2010|for reaction 13 Stenrup, Larson & Elander|2009|for reaction 



IS 



and 

%3]/^T3~~0 03 to within~50% for all temperatures 100 < f < 10 4 K. If we com- 
pare this value with the residual f ractional ionization of the intergalactic medium 
(IGM) at this epoch, x — 2 x 10~ 4 (|Schleicher et al. 



2008 ), we see that mutual neu- 



tralization is unimportant within the very first star-forming minihalos. It becomes 
important once larger minihalos, with virial temperatures r v ; r ~ 10 4 K or above, be- 
gin to form, as in these minihalos, substantial collisional ionization of the gas can 
occur, leading to an initial fractional ionization much higher than the residual value 
in the IGM. It also becomes an important process within the "fossil" HII regions left 



behind by the first generation of massive stars (|Oh & Haiman| |2003 Nagakura & 



Omukai 2005 Glover, Savin & Jappsen 2006; Kreckel et al. , 2010). 



1 A group lead by X. Urbain at the Universite Catholique de Louvain has recently made new 
experimental measurements of the rate of this reaction at low temperatures, but at the time of 
writing, the results of this work remain unpublished 
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Although H2 is by far the most abundant primordial molecule, it is actually not a 
particularly efficient coolant. The H2 molecule has no dipole moment, and so dipole 
transitions between its excited rotational and vibrational levels are forbidden. Al- 
though radiative transitions between levels do occur, they are quadrupole transitions 
and the associated transition rates are small. In addition, application of the Pauli 
exclusion principle to the hydrogen molecule shows that it must have two distinct 
states, distinguished by the nuclear spin of the two hydrogen nuclei: para-hydrogen, 
in which the nuclear spins are parallel, and which must have an even value for the 
rotational quantum number 7, and ortho-hydrogen, which has anti-parallel nuclear 
spins and an odd value for 7. Radiative transitions between ortho-hydrogen and 
para-hydrogen involve a change in orientation of the spin of one of the nuclei and 
are therefore strongly forbidden. As a result, the least energetic rotational transition 
of H2 that has any significant probability of occurring is the transition between the 
7 = 2 and 7 = rotational levels in the vibrational ground-state of para-hydrogen. 
This transition has an associated energy £20/^ — 512 K. The H2 molecule there- 
fore has large energy separations between the ground state and any of the accessible 
excited rotational or vibrational state^] and has only weak radiative transitions be- 
tween these states. 

These features of the H2 molecule have two important consequences. First, it 
becomes a very inefficient coolant at temperatures T -C E20/K, as it becomes al- 
most impossible to collisionally populate any of the excited states. The minimum 
temperature that can be reached solely with H2 cooling depends somewhat on the 
H2 abundance and the time available for cooling, but typically 7 m j n ~ 150-200 K. 
Second, its rotational and vibrational levels reach their local thermodynamic equi- 
librium (LTE) level populations at a relatively low density, n cl ; t ~ 10 4 cm~ 3 . This 
means that at densities n ^> n cr ; t , the H2 cooling rate scales only linearly with den- 
sity and the cooling time due to H2 becomes independent of density. Since other 
important timescales, such as the free-fall collapse time of the gas, continue to de- 
crease with increasing density, the implication is that H2 becomes an increasingly 
ineffective coolant as one moves to higher densities. 

For these reasons, primordial molecules or molecular ions that do not share these 



drawbacks have attracted a certain amount of attention. In an early study, Lepp & 
|Shull| ( |1984| l suggested that deuterated hydrogen, HD, and lithium hydride, LiH, 
may both be significant coolants in primordial gas. More recently, work by Yoshida 



et al. ( 2007 1 h as suggested tha t H, m ay be an important coolant in some circum- 
stances, while Glover & Savin (2006) show that H^ is also worthy of attention. In 



practice, the only one of these molecules or ions that has proved to be important is 
HD. Detailed modelling of the chemistry of hthium in primordial gas (e.g. |Stancil,| 
|Lepp & Dalgarno 1996 1 Mizusawa, Omukai & Nishi 2005 ) has shown that LiH is 



efficiently destroyed by the reaction 

LiH + H^Li + H 2 , (20) 



2 For comparison, note that the energy separation between the / = and / = 1 rotational levels of 
CO is roughly 5 K. 
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and that only a small fraction of the available lithium (which itself has an a bundance 
of only 5 x 10~ 10 relative to hydrogen; see Cyburt, Fields & Olive 2008 ) is ever 
incorporated into LiH. Cooling from the molecular ion Hj was re-examined by 
Glover & Savin | (12009 ), who showed that the collisional excitation rates cited by 
Galli & Palla (1998 ) and used as a basis for the fits given in Yoshida et al. (2007) 
were a factor of ten too large, and that if the correct rates are used, Hj cooling 
is no longer important. Finally, Glover & Savin (2009) also examined the possible 
role played by H^ cooling in considerable detail, but found that even if one makes 
optimistic assumptions regarding its formation rate and collisional excitation rate, it 
still contributes to the total cooling rate at the level of only a few percent, and hence 
at best is a minor correction term. 

These studies leave HD as the only viable alternative to H2 as a coolant of pri- 
mordial gas. HD has a small, but non-zero dipole moment, giving it radiative tran- 
sition rates that are somewhat larger than those of H2, resulting in a critical density 
n crit ~ 10 6 cm~ 3 . Unlike H2, it is not separated into ortho and para states, and so the 
lowest energy transition accessible from the ground state is the 7=1 to J — ro- 
tational transition, with an energy Eio/k = 128 K. Although the c osmological ratio 



of deuterium to hydrogen is small [D/H = (2.49 ±0.17) x lO MCyburt, Fields & 
Olive|2008) , the ratio of HD to H2 can be significantly boosted in low temperature 
gas by chemical fractionation. The reaction 



H 2 + D+ ^HD + H^ 



(21) 



that converts H2 into HD is exothermic and so proceeds rapidly at all temperatures, 
while the inverse reaction 



HD + H+->H 2 +D^ 



(22) 



is endothermic and so proceeds very slowly at low temperatures. In equilibrium, 
these two reactions produce an HD-to-H2 ratio given by 



™=2exp(^)[D/H], 
*H, \ T J 



(23) 



where [D/H] is the cosmological D:H ratio. Together, these factors render HD a 
much more effective coolant than H2 in low temperature gas. 

In practice, for HD cooling to take over from H2 cooling, the gas must already be 
fairly cold, with T ~ 150 K (Glover 2008) >, and temperatures this low are typically 
not reached during the collapse of the first star-forming minihalos, meaning that 



HD remains a minor coolant (Bromm, Coppi & Larson 2002 1. However, there are a 
number of situations, typically involving gas with an enhanced fractional ionization 
in which HD cooling does become significant (see e.g. Nakamura & Umemura 

2006 



2002 Nagakura & Omukai 2005 1 Johnson & Bromm 



McGreer & Bryan||2008||Greif et al.||2008||Kreckel et al.||20lO) T 



lYoshid a etaL] |2007 



>7; 
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1.3 The minimum mass scale for collapse 



The relative simplicity of the chemistry discussed in the previous section allows one 
to construct a very simple model that captures the main features of the evolution of 
the H2 fraction within low density gas falling into a dark matter minihalo. We start 
by assuming that radiative recombination is the only process affecting the electron 
abundance, and writing the rate of change of the electron number density as 

d« e 

-r- = -&rec«e«H+) (24) 
at 

where n e is the number density of electrons, n H + is the number density of protons, 
and k lec is the recombination coefficient. If we assume that ionized hydrogen is 
the only source of free electrons, implying that n e = «h+> an d tnat tne temperature 
remains roughly constant during the evolution of the gas, then we can solve for the 
time evolution of the electron fraction: 

(25) 



1 + k mc ntxo ' 



where x = n e /n, n is the number density of hydrogen nuclei, and xq is the initial 
value of x. We next assume that all of the H2 forms via the H pathway, and that 



mutual neutralization of H with H + (reaction 19 1 is the only process competing 



with associative detachment (reaction 13 1 for the H ions. In this case, we can write 
the time evolution of the H2 fraction, xn 2 = «H 2 / w > as 



d*H 2 
dt 



^npmpAD, (26) 



where ^jjjis the rate coefficient of reaction 12 the formation of H by the radiative 
association of H and e~, and pad is the probability that any given H ion will be 
destroyed by associative detachment rather than by mutual neutralization. Given our 
assumptions above, this can be written as 

tern , 07 n 

Pad = ^fr. , (27) 

where Aro] is the rate coefficient for reaction [13] and Arrg] is the rate coefficient for 
reaction [19] If we again assume that n e = « H +, and in addition assume that «h — «. 
then the expression for pad can be simplified to 

PAD=(l + g*)" (28) 



Substituting this into Equation 26 we obtain 



*** -HnmU + tex) ■ (29) 



^ V tea 
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If the initial fractional ionization xq <C frnl/ ^ fT9t then the term in parentheses is of 
order unity and this equation has the approximate solution 

xn 2 = ^2 In (1 + k mc nx t ) , (30) 

= ^ln(l+t/ftec), (31) 

where f rec = 1/ (£rec«*o) is the recombination time. The growth of the H2 fraction 
is therefore logarithmic in time, with most of the H2 forming within the first few 
recombination times. In the more complicated case in which xq is comparable to or 
larger than /fy}]/lfi9[ but still significantly less than unity (so that «h ~ i)> the H2 
fraction is given instead by 

^fla, ( 1 +- y o A Eii/tei+ f A" 



x h 2 = T^ln ( ) ■ (32) 



From this analysis, we see that the main factor determining the final H2 abun- 
dance is the ratio /^i2\/k mc , since for times of the order of a few recombination times, 



the logarithmic term in Equation 32 is of order unity, implying that the final H2 abun- 
dance is at most a factor of a few times tel /fcrec- ^ we use tne simple power-law fits 
to Afngand k lec given by |Hutchins| ( |l976[ >, namely 1.83 x 1(T 18 r°- 8779 cm 3 s" 1 



and k lec = 1.88 x 10 lu T u b44 cm 3 s , then we can write the ratio of the two rate 
coefficients as 

^ 10 -8 r 1.5219 (33) 

The amount of H2 produced is a strong function of temperature, but is of the order 
of a few times 10~ 3 for temperatures of a few thousand Kelvin. We see therefore 
that the formation of H2 via H~ never results in a gas dominated by H2, as the H2 
abundance always remains much smaller than the abundance of atomic hydrogen. 

Given this simple model for the amount of H2 that will form in the gas, the 
obvious next step is to compare this to the amount of H2 that is required to cool the 
gas efficiently. In order to determine the H2 fraction necessary to significantly cool 
gas with a temperature T within some specified fraction of the Hubble time - say 
20% of ?h - we can simply equate the two timescales, and solve for the H2 fraction. 
Using our previous definition of the cooling time, we have 

1 nMkT :0.2r H , (34) 



J- 1 A (7>h 2 



where we have assumed that H2 is the dominant coolant and have written the cooling 
rate per unit volume in terms of Aq, the cooling rate per H2 molecule. Rearranging 
this equation, using the fact that when the H2 fraction and the ionization level are 
low, 7=5/3 and n tot = ( 1 + 4xn e )», where ^He is the fractional abundance of helium 
(given by xn e = 0.083 for primordial gas), we obtain 
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XH 2! req = 1.38xlO- 15 ^— (35) 

In the high-redshift limit where ?h — Hq 1 £1,„ ^ 2 (1 + z)~ 3 / 2 , this becomes 



, H2 , req = 5.2xlCT-—^—) , (36) 



where we have used values for the cosmological parameters taken from Komatsu 



et al. ( 201 1) . Collisions of H2 with a number of different species contribute to Aq, 



as explored in Glo ver & Abel| ( |2"008| >, but in the earliest minihalos, the dominant 
contributions come from collisions with H and He. Aq is therefore given to a good 
approximation by 

Aq = A H n H + A He «He- (37) 

Simple fits for the values of Ah and An e as a function of temperature can be found 
in |Glover & Abel| | |2008l >. 

An illustration of the likely size of XH 2 , re q is given in Figure [2] In this Figure, 
we plot xn,,req as a function of temperature, evaluated for three different redshifts: 
z = 20, 30 and 40. In computing these values, we have assumed that the mean density 
of the gas in the minihalo is given by p = A c p&o> where po is the cosmological 
background density of baryons. In the Figure, we also show the actual H2 fraction 
produced in the gas, XH 2 , ac t, as a function of temperature at times equal to 1, 5 and 
10 recombination times, and where we have taken Xq -C ^nj/^Jij} 

Figure|2]demonstrates that the amount of H2 produced in the gas is a strongly in- 
creasing function of temperature, while the amount required to bring about efficient 
cooling of the gas is a strongly decreasing function of temperature. This means that 
for any given choice of comparison time t and redshift z, we can identify a critical 
temperature r cr i t , such that gas with T > T a i t will cool within a small fraction of a 
Hubble time, while gas with T < r cl ; t will not. Moreover, because XH 2 , ac t an d XH 2 ,req 
are both steep functions of temperature, but are relatively insensitive to changes in 
t or z, the value of T CI it that we obtain is also relatively insensitive to our choices 
for t or z. We find that r cr it ~ 1000 K, and that at this temperature, the H2 frac- 
tion required t o provide efficient co oling lies somewhere between a few times 10~ 4 



and 10 3 (c.f.|Tegmark et al. 1997 who come to a similar conclusion using a very 



similar argument). If we convert this critical virial temperature into a corresponding 



critical minihalo mass using Equation 10 we find that 



M^-exloV 1 ^) 3/2 .Q,„ 1/2 (^) ' M . (38) 

This mass scale is illustrated by the dotted line in Figure [T] At high redshift, it is 
smaller than the filter mass scale corresponding to Vb c = C7 v bc> demonstrating that 
at these redshifts, it is the streaming of the gas with respect to the dark matter that 
is the main process limiting the formation of Population III stars. Below a redshift 
of around 40, however, M cr ; t becomes the larger mass scale, implying that at these 
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Fig. 2 Comparison of the fractional abundance of H2 produced with our simple toy model for the 
chemistry (solid lines) versus the quantity of H2 required in order to cool the gas within 20% of 
a Hubble time (dashed lines). From bottom to top, the solid lines correspond to the H2 fraction 
produced at times t = 1,5 and 10f rec , respectively, where ? rec is the recombination time, and the 
dashed lines correspond to the H2 fraction required at redshifts z = 40, 30 and 20, respectively. 
We see that the minimum temperature that the gas must have in order to be able to cool within a 
fraction of a Hubble time - indicated by the point at which the lines cross - is relatively insensitive 
to our choices for t and z. 



lower redshifts, there will be a population of small minihalos that contain a signif- 
icant gas fraction, but that do not form stars, because their gas is unable to cool in 
less than a Hubble time. These small starless minihalos may be important sinks for 
ionizing photons during the epoch of reionization (Haiman, Abel & Madau 2001 ). 



To conclude our discussion of the first star-forming minihalos, we should men- 
tion one potentially important effect not taken into account in the analysis above. 
This is the influence of ongoing minor mergers and accretion on the thermal balance 
of the gas. Although major mergers occur only once per Hubble time, on average, 
minor mergers occur far more frequently, and act to stir up the gas, thereby heating 
it and lengthening the time required for it to cool. This phenomenon was noted by 
|Yoshida et al.| ( [2003] l in their cosmological simulations of the formation of the first 
star-forming minihalos. Yoshida et al.| ( 2003] l show that in spite of the approxima- 
tions made in its derivation, Equation [38] gives a reasonable guide to the minimum 
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mass of the minihalos that contain gas that can cool effectively. However, they also 
find that there are some minihalos with M > M cr ; t in which the gas does not cool. 
They show that these minihalos have higher mass accretion rates than minihalos 
of the same mass in which cooling does occur, and hence ascribe the suppression 
of cooling to the effects of dynamical heating by the ongoing accretion and minor 
mergers. This effect was also treated more recently by |Wang & Ab el (2008 ), who 
show that it can be included into the simple thermal model described above by the 
addition of a heating term describing the effects of mergers and accretion. They 
show that if one writes this heating term as 



k dT, 



vir 



y- 1 df 



(39) 



then one can relate the rate of change of the virial temperature to the mass growth 
rate of the minihalo in a relatively simple fashion. 



2 Gravitational collapse and the formation of the first protostar 

As the analysis in the previous section has shown, gas in minihalos with virial tem- 
peratures greater than about 1000 K (corresponding to masses M ~ 1O 6 M ) can 
form enough H2 to cool within a small fraction of a Hubble time. This reduces 
the pressure and allows the gas to collapse further under the influence of its own 
self-gravity. As it does so, the value of the Jeans mass decreases. Many early stud- 
ies of the formation of primordial stars assumed that as the Jeans mass decreases 
and the gas becomes more gravitationally unstable, it begins to undergo hierarchical 
gravitational fragmentation in a manner similar to that envisaged by |Hoyle] ( |1953[ ), 
with the result that at any given moment, the mean fragment mass is approximately 
equal to the local Jeans mass (see |Glover 2005 1 for a historical summary of these 



models). In this picture, one could predict the final mass of the first stars simply by 
studying the evolution of the Jeans mass. Moreover, since the minimum Jeans mass 
reached during the collapse can be estimated with reasonable accuracy on purely 
thermodynamic al grounds ( |Rees| [T9"76| |Low & Lynd en-Bell| [T976 ), in this view of 
Population III star formation, the dynamics of the gas is of secondary importance. 
Around ten years ago, however, it first became possible to model the coupled chem- 
ical, dynamical and thermal evolution of the gas within a primordial minihalo using 



high resolution 3D numerical simulations (Abel, Bryan & Norman 2000 2002 
Bromm, Copp i & Larson| [T999 2002). These studies showed that the picture out- 
lined above is wrong: the gas does not undergo hierarchical fragmentation, and so 
one cannot predict the masses of the first stars simply by studying the evolution 
of the Jeans mass. These high resolution numerical simulations, and the many that 



have followed them (e.g. Yoshida et al.| |2006[ |Q'Shea & Normanj |2007[ |McGreer 



|& Bryan[|2008| to name but a few), have for the first time given us a clear picture of 
exactly how gravitational collapse proceeds within one of these early minihalos. In 
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the next section, we will discuss the sequence of events that occur as we follow the 
collapse from the minihalo scale all the way down to the scale of a single Popula- 
tion III protostar. Following that, in Sections 2.2 and 2.3 we discuss two of the main 
uncertainties remaining in our model for the formation of the first Pop. Ill protostar: 
the role played by heating and ionization arising from dark matter self-annihilation 
(Section|2~2|) and the role played by magnetic fields (Section|2~3|). 



2.1 Thermal and chemical evolution of the gas during collapse 
2.1.1 Initial collapse 

As gas falls into the minihalo from the intergalactic medium, it is shock-heated to a 
temperature close to r v ; r . In the post-shock gas, the electron fraction decreases due 



to radiative recombination, but at the same time H2 forms, primarily via reactions 12 
and[T3] As we have already seen, the H2 fraction evolves logarithmically with time, 
with most of the H2 forming within the first few recombination times. As the H2 
fraction increases, so does its ability to cool the gas, and so the gas temperature 
slowly decreases, reducing the pressure and allowing the gas to collapse to the centre 
of the minihalo. 

At this point, the evolution of the gas depends upon how much H2 it has formed. 
There are two main outcomes, and which one occurs within a given minihalo de- 
pends primarily on the initial ionization state of the gas. 



The low ionization case 

During the formation of the very first Population III stars (also known as Population 
III. 1, to use the terminology introduced by Tan & McKee 2008), the initial frac- 



tional ionization of the gas is the same as the residual ionization in the intergalactic 
medium, i.e. xo ~ 2 x 10~ 4 . In this case, the amount of H2 that forms in the gas is 
typically enough to cool it to a temperature of T ~ 200 K but not below. At this tem- 
perature, chemical fractionation has already increased the HD/H2 ratio by a factor 
of 20 compared to the cosmic deuterium-to-hydrogen ratio, and as a consequence, 
HD is starting to become an important coolant. However, the amount of HD that 



forms in the gas is not enough to cool it significantly below 200 K (Bromm, Coppi 



|& Larson] |2002[ ), and H2 continues to dominate the cooling and control the further 
evolution of the gas. In this scenario, the collapse of the gas is greatly slowed once 
its temperature reaches 200 K and its density reaches a value of around 10 4 cm -3 , 
corresponding to the critical density n cr j t , at which the rotational and vibrational 
level populations of H2 approach their local thermodynamic equilibrium (LTE) val- 
ues. At densities higher than this critical density, the H2 cooling rate per unit volume 
scales only linearly with n (compared to a quadratic dependence, Ar 2 « n 2 at lower 
densities), while processes such as compressional heating continue to increase more 
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rapidly with n. As a result, the gas temperature begins to increase once the density 
exceeds w cr j t - 

Gas reaching this point in the collapse enters what Bromm, Coppi & Larson 



(2002) term a "loitering" phase, during which cold gas accumulates in the centre 
of the halo but only slowly increases its density. This loitering phase ends once the 
mass of cold gas that has accumulated exceeds the local value of the Bonnor-Ebert 
mass (Bonnor 1956; Ebert 1955), given in this case by (Abel, Bryan & Norman, 

[20021 

M BE ^ 40T 3 / 2 n~ 1/2 M®, (40) 

which for n ~ 10 4 cm~ 3 and T ~ 200 K yields M BE ~ lOOOM 0Once the mass of 
cold gas exceeds Mbe, its collapse speeds up again, and becomes largely decoupled 
from the larger-scale behaviour of the gas. The next notable event to occur in the 
gas is the onset of three-body H2 formation, which is discussed in the next section. 



The high ionization case 

If the initial fractional ionization of the gas is significantly higher than the residual 
fraction in the IGM, then a slightly different chain of events can occur. A larger 
initial fractional ionization implies a shorter recombination time, and hence a log- 
arithmic increase in the amount of H2 formed after a given physical time. An in- 
crease in the H2 fraction allows the gas to cool to a slightly lower temperature, and 
hence boosts the HD abundance in two ways: the lower temperature increases the 
HD/H2 ratio produced by fractionation, and the H2 fraction itself is larger, so any 
given HD/H2 ratio corresponds to a higher HD abundance than in the low ionization 
case. If the addition ionization allows enough H2 to be produced to cool the gas 
to T ~ 150 K (which requires roughly a factor of three more H2 than is required 
to reach 200 K), then chemical fractionation increases the HD abundance to such 
an extent that it takes over as the dominant coolant ( |Glover| |2008| ). This allows 
the gas to cool further, in some cases reaching a temperature as low as the CMB 
temperature, 7cmb (e.g. Nakamura & Umemura, 2002; Nagakura & Omukai, 2005 



Johnson & Bromm 2006 Yoshida et al. 2007) |McGreer & Bryan| |2008| |Kreckel 



5 
el 



et al.| 20 1 Oj l . The higher critical density of HD, n c rit,HD ~ 10 b cm means that the 



gas does not reach the loitering phase until much later in its collapse. Once the gas 
does reach this phase, however, its subsequent evolution is very similar to that in the 
low-ionization case discussed above. Cold gas accumulates at n ~ n cr ; t until its mass 
exceeds the Bonnor-Ebert mass, which in this case is Mbe ~ 40 Mq if T = 100 K 
and n = 10 6 cm -3 . Once the gas mass exceeds Mbe, the collapse speeds up again, 
and the gas begins to heat up. Aside from the substantial difference in the size of 
Mbe, the main difference between the evolution of the gas in this case and in the low 



3 Discussions of Population III star formation often refer to the cold clump of gas at the centre of 
the minihalo as a "fragment", and speak of Mbe as the "fragmentation mass scale", but in the case 
of the very first generation of star-forming minihalos, this is actually something of a misnomer, as 
very seldom does more than one "fragment" form in a given minihalo 
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ionization case lies in the fact that in the high ionization case, the gas reheats from 
T ~ 200 K or below to T ~ 1000 K much more rapidly than in the low ionization 
case. As we shall see later, this period of rapid heating has a profound influence on 
the ability of the gas to fragment. 

Several different scenarios have been identified that lead to an enhanced frac- 
tional ionization in the gas, and that potentially allow the gas to reach the HD- 
dominated regime. Gas within minihalos with T v - lr > 9000 K will become hot enough 
for collisional ionization of hydrogen to supply the necessary electrons. However, as 



halos of this size will typically have at least one star-forming progenitor (Johnson, 
2008 ), it is questionable whether many Pop. Ill stars will form in 



Greif & Bromm 



such minihalos, as we would expect the gas in most of them to have been enriched 
with metals by one or more previous episodes of star formation. 

Another possibility that has attracted significant attention involves the gas in the 
minihalo being drawn from a "fossil" HII region, i.e. a region that was formerly 
ionized by a previous Population III protostar but has now recombined (see e.g. Oh 
& Haimanj |2lX)3| |Nagakura & Omukaij [2005] |Yoshida et al) [2007] ). Many studies 



have shown that the volume of the IGM ionized by a single massive Pop. Ill star is 
significantly larger than the volume that is enriched by the metals produced in the 
supernova occurring at the end of the massive star's life (see e.g. the recent treatment 
by Greif et al. 2010 or Ciardi & Ferrara ( 2005 \ for a summary of earlier work). It is 



therefore possible that a significant number of Population III stars may form in such 
conditions. 

A final possibility is that the required ionization can be produced by a flux of X- 
rays or high energy cosmic rays. Although X-ray ionization was initially favoured 
as a means of raising the ionization level of the gas, and hence promoting H2 forma- 
tion ([Haiman, Abel & Re es , 2000), more recent work has shown that if one consid- 
ers realistic models for the X-ray background that also account for the simultaneous 
growth of the soft UV background, then one finds that UV photodissociation of 
the H2 is a more important effect, and hence that the growth of the radiation back- 
grounds almost always leads to an overall reduction in the amount of H2 produced 
(Glover & Brand 2003 Machacek, Bryan & Abel| |2Q03]>. Cosmic ray ionization 



JM2_ 

may therefore prove to be the more important effect (Stacy & Bromm] |2007[|Jasche 



|Ciardi & Ensslin| [2007 ), although we still know very little about the likely size of 
the cosmic ray ionization rate in high redshift minihalos. 



2.1.2 Three-body H2 formation 



Once the collapsing gas reaches a density of around 10 —10 cm , its chemical 
makeup starts to change significantly. The reason for this is that at these densities, 



the formation of H2 via the three-body reactions (Palla, Salpeter, & Stahler 1983 ) 



H + H + H H 2 +H, 
H + H + H 2 ^H 2 +H 2 , 
H + H + He -» H 2 +He, 



(41) 
(42) 
(43) 
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starts to become significant. These reactions quickly convert most of the hydrogen 
in the gas into H2. At the same time, however, they generate a substantial amount 
of thermal energy: every time an H2 molecule forms via one of these three-body 
reactions, its binding energy of 4.48 eV is converted into heat. A simple estimate of 
the relative sizes of the compressional heating rate and the three-body H2 formation 
heating rate helps to demonstrate the importance of the latter during this stage of 
the collapse. Let us consider gas at a density n = 10 8 cm" 3 that has a temperature 
T = 1000 K, collapsing at a rate such that dn/dt =«/%, where % is the gravitational 
free-fall time. In these conditions, the compressional heating rate is given by 

A pdv ~ 1.25 x 10" 3 V /2 7\ (44) 
= 1.25 x 10" 24 ergs" 1 cm" 3 , (45) 

while the three-body H2 formation heating rate has the value 

A 3b ^ 3.9 x 10" 40 r" 1 n 3 x 3 1 , (46) 
= 3.9 x 10" 19 4 ergs" 1 cm" 3 , (47) 



where we have adopted the rate coefficient for reaction 41 given in Palla, Salpeter, 



& Stahler ( 1983[ ). Comparing the two heating rates, we see that three-body H2 for- 



mation heating dominates unless jch is very small (i.e. unless the gas is almost fully 
molecular). Therefore, even though the abundance of H2, the dominant coolant dur- 
ing this phase of the collapse, increases by more than two orders of magnitude, the 
gas typically does not cool significantly, owing to the influence of this three-body 
H2 formation heating. Indeed, the temperature often actually increases. 

One major uncertainty that remains in current treatments of this phase of the 
collapse of the gas is exactly how quickly the gas becomes molecular. Although re- 
action|4T]is the dominant source of H2 at these densities, the rate coefficient for this 
reaction is poorly known, with published values differing by almost two orders of 



magnitude at 1000 K, and by an even larger factor at lower temperatures (Glover 



2008 Turk et al. |2011 1. The effects of this uncertainty have recently been studied 
by |Turk et al. ( 201 l| l. They show that it has little effect on the density profile of 



the gas, and only a limited effect on the temperature profile. However, it has much 
more significant effects on the morphology of the gas and on its velocity structure. 
Simulations in which a high value was used for the three-body rate coefficient show 
find that gas occurs more rapidly, and that the molecular gas develops a much more 
flattened, filamentary structure. Significant differences are also apparent in the infall 
velocities and the degree of rotational support. Turk et al. ( 2011| l halt their simula- 



tions at the point at which a protostar first forms, and so do not directly address the 
issue of whether these differences continue to have an influence during the accre- 



tion phase, and whether the affect the ability of the gas to fragment (see Section 3.2 
below). A follow-up study that focussed on these issues would be informative. 
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2.1.3 Optically-thick line cooling 



The next important event occurs at a density of around 10 10 cm i , when the main 



rotational and vibrational lines of H2 start to become optically thick ( Ripamonti & 



Abel 2004). The effect of this is to reduce the efficiency of H2 cooling, leading to a 
continued rise in the gas temperature. In one-dimensional simulations (e.g. Omukai 
& Nishi| |1998| |Omukai et al] |1998[ |Ripamonti et al.[ |2002| >, it is possible to treat 
optically thick H2 cooling accurately by solving the full radiative transfer problem. 
These models show that although the optical depth of the gas becomes large at 
frequencies corresponding to the centers of the main H2 emission lines, the low 
continuum opacity of the gas allows photons to continue to escape through the wings 
of the lines, with the result that the H2 cooling rate is suppressed far less rapidly as 



the collapse proceeds than one might at first expect (see Omukai et al. 1998 for a 
detailed discussion of this point). 

In three-dimensional simulations, solution of the full radiative transfer problem is 
not currently possible, due to the high computational expense, which has motivated 
a search for simpler approximations. There are two such approximations in current 
use in simulations of Population III star formation. The first of these was introduced 



by Ripamonti & Abel ( 2004 1. They proposed that the ratio of the optically thick and 



optically thin H2 cooling rates, 



A 



H 2 . thick 
Ah,, thin ' 



could be represented as a simple function of density: 



fx 



-0.45 



(48) 



(49) 



where no = 8 x 10 cm . They showed that this simple expression was a good ap- 



proximation to the results of the full radiative transfer model used by Ripamonti 
|et al.| ( |2002[ ), and suggested that this approximation would be useful for extending 
the results of three-dimensional simulations into the optically thick regime. How- 
ever, they also noted that it may only be accurate while the collapse remains ap- 
proximately spherical, as the one-dimensional model on which it is based assumes 
spherical infall. 

An alternative approach was introduced by Yosh ida et al.| (|2006). They compute 
escape probabilities for each rotational and vibrational line of H2 using the standard 
Sobolev approximation (Sobolev, 1960). In this approximation, the optical depth at 
line centre of a transition from an upper level u to a lower level I is written as 



(50) 



where a u \ is the line absorption coefficient and L s is the Sobolev length. The ab- 
sorption coefficient a„i can be written as 
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Oui = —. — niBh, 



<Kv„/), 



(51) 



where E u i is the energy difference between the two levels, n\ is the number density 
of H2 molecules in the lower levels, B\ u is the usual Einstein B coefficient, and 
(v u i) is the line profile at the centre of the line. The Sobolev length is given by 



i'th 



|dv r /dr| 



(52) 



where v t h is the thermal velocity of the H2 and |dv r /dr| is the size of the velocity 
gradient along a given line of sight from the fluid element of interest. Given T H /, the 
escape probability for photons emitted in that direction then follows as 



1 -exp(-T H/ ) 

*ul 



(53) 



To account for the fact that the velocity gradient may differ along different lines of 
sight from any particular fluid element, |Yoshida et al. ( 2006| l utilize a mean escape 
probability given by 

^ = MMft, (54 ) 

where /5 X , j5 y and j5 z are the escape probabilities along lines of sight in the x, y and 
Z directions, respectively. Finally, once the escape probabilities for each transition 
have been calculated, the optically thick H2 cooling rate can be computed from 



^H 2 .thick — ^E u iPuiA u in u , 

u.l 



(55) 



where A„/ is the Einstein A coefficient for the transition from u to I and n u is the 
population of the upper level u. 

Strictly speaking, the Sobolev approximation is valid only for flows in which the 
Sobolev length L s is much smaller than the characteristic length scales associated 
with changes in the density, temperature or chemical makeup of the gas, a require- 
ment which is easy to satisfy when the velocity gradient is very large, but which is 
harder to justify in the case of Population III star formation, since the collapse speed 



is typically comparable to the sound-speed. Nevertheless, Yoshida et al. (2006) show 
that the optically thick H2 cooling rates predicted by the Sobolev approximation 
are in very good agreement with those computed in the one-dimensional study of 
Omukai & Nishi ( 1998) by solution of the full radiative transfer problem. 

Little work has been done on comparing these two approaches to treating 



optically-thick H2 cooling. This issue was addressed briefly in Turk et al. (2011 1 



who showed that the two approximations yielded similar values for f T for densities 



n < 10 cm during the initial collapse of the gas, with differences of at most a 
factor of two. However, as yet no study has examined whether this good agreement 
persists past the point at which the first protostar forms. 
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2.1.4 Collision-induced emission 



A further significant point in the collapse of the gas is reached once the number 
density increases to n~ 10 14 cm~ 3 . At this density, a process known as collision- 
induced emission becomes important. Although an isolated H2 molecule has no 
dipole moment, and can only emit or absorb radiation through quadrupole transi- 
tions, when two H2 molecules collid^jthey briefly act as a kind of "supermolecule" 
with a non-zero dipole moment for the duration of the collision. This supermolecule 
can therefore absorb or emit radiation through dipole transitions, which have much 
higher transition probabilities than the quadrupole transitions available to isolated 
H2. If radiation is absorbed, this process is termed collision-induced absorption; if 
it is emitted, then we refer to the process as collision-induced emission (CIE). A 
detailed discussion of the phenomenon can be found in Frommhold ( 1993 ). 

Collision-induced emission can in principle occur in gas of any density, but the 
probability of a photon being emitted in any given collision is very small, owing to 
the short lifetime of the collision state (At < 10~ 12 s at the temperatures relevant for 
Pop. Ill star formation; see Ripamonti & Abel 2004). For this reason, CIE becomes 
an important process only at very high gas densities. Another consequence of the 
short lifetime of the collision state is that the individual lines associated with the 
dipole transitions become so broadened that they actually merge into a continuum. 
This is important, as it means that the high opacity of the gas in the rovibrational 
lines of H2 does not significantly reduce the amount of energy that can be radiated 
away by CIE. Therefore, once the gas reaches a sufficiently high density, CIE be- 



comes the dominant form of cooling, as pointed out by several authors ( |Omukai & 
MsMl[T998l|Ripamonti etaL] [2002) [Ripamonti & Abel|[2004| . 

The most detailed study of the effects of CIE cooling on the collapse of primor- 



dial gas was carried out by Ripamonti & Abel (2004). They showed that CIE cooling 
could actually become strong enough to trigger a thermal instability, However, the 
growth rate of this instability is longer than the gravitational free-fall time, meaning 
that it is unlikely that this process can drive fragmentation during the initial collapse 
of the gas. 



2.1.5 Cooling due to H 2 dissociation 



The phase of the collapse dominated by CIE cooling lasts for only a relatively short 
period of time. The gas becomes optically thick in the continuum once it re aches a 



density n ~ 10 16 cm~ J (|Omukai & Nishrj |1998[|Ripamonti & Abel| |2004|), which 



strongly suppresses any further radiative cooling. Once this occurs, the gas temper- 
ature rises until it reaches a point at which the H2 begins to dissociate. At these 
densities, this occurs at a temperature T ~ 3000 K. Once this point is reached, the 
temperature rise slows, as most of the energy released during the collapse goes into 
dissociating the H2 rather than raising the temperature. As it takes 4.48 e V of energy 



4 A similar process can also occur during collisions of atomic hydrogen or atomic helium with 
H2, but it is the H2-H2 case that is the most relevant here 
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to destroy each H2 molecule, this H2 dissociation phase continues for a while. How- 
ever, it comes to an end once almost all of the H2 has been destroyed, at which point 
the temperature of the gas begins to climb steeply. The thermal pressure in the inte- 
rior of the collapsing core rises rapidly and eventually becomes strong enough to halt 
the collapse. At the point at which this occurs, the size of the dense core is around 
0.1 AU, its mass is around 0.01 M and its mean density is of order 10 20 cm -3 
( Yoshida, Omukai & Hernquist 2008[ ). It is bounded by a strong accretion shock. 
This pressure-supported, shock-bounded core is a Population III protostar, and its 
later evolution is discussed in Section[3]below. 



2.2 Dark matter annihilation 

One complication not accounted for in the models of Pop. Ill star formation de- 
scribed above is the role that may be played by dark matter annihilation. The nature 
of dark matter is not yet understood, but one plausible candidate is a weakly inter- 
acting massive particle (WIMP) - specifically, the lightest supersymmetric particle 
predicted in models based on supersymmetry. The simplest supersymmetry models 
predict that this WIMP has an annihilation cross-section (crv) ~3x 10~ 26 cm 2 , a 
mass within the range of 50 GeV to 2 TeV, and a cosmological density consistent 
with the inferred density of dark matter ( |Spolyar et al.| |2008| l. The rate per unit 
volume at which energy is produced by dark matter annihilation can be written as 
Gann = ( c7y )pl / m x, where px is the mass density of dark matter and nix is the mass 
of a single dark matter particle. For a plausible particle mass of 100 GeV, and a dark 
matter density equal to the cosmological background density of dark matter, this ex- 
pression yields a tiny heating rate, Q atm ~6x 10~ 62 (1 +z) 6 i2, 2 1 /i 4 ergcm~ 3 s _1 , even 
before one accounts for the fact that much of the annihilation energy is released in 
the form of energetic neutrinos or gamma-rays that couple only very weakly with the 
intergalactic gas. WIMP annihilation therefore plays no significant role in the evo- 
lution of the intergalacti c medium while the WIMPs remain uniformly distributed 



(Myers & Nusser 2008). However, the p x density dependence of the heating rate 
means that it can potentially become significant in regions where the dark matter 
density is very high. 



Spolyar et al. (2008 ) proposed that one situation in which the heating from dark 



matter annihilation could become important would occur during the formation of 
the very first Population III protostars. They assumed that any given star-forming 
minihalo would form only a single Pop. Ill protostar, and that this protostar would 
form at the center of the minihalo. As the gas collapsed at the center of the minihalo, 
its increasing gravitational influence would bring about a local enhancement of the 
dark matter density, via a process known as adiabatic contraction. The basic idea 
underlying this is very simple. For a collisionless particle on a periodic orbit, the 
quantity § pdq, where p is the conjugate momentum of coordinate q, is an adiabatic 
invariant, i.e. a quantity that does not vary when the gravitational potential varies, 
provided that the rate of change of the potential is sufficiently slow. If p represents 
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the angular momentum of a particle on a circular orbit of radius r within some 
spherically symmetric mass distribution, then one can show that the quantity rM(r) 
is constant for that particle, where M(r) is the mass enclosed within r, so long as this 
enclosed mass changes on a timescale that is long compared to the orbital period. 
Spolyar et al. ( 2008|l show that if o ne starts with a simple NFW profile for the dark 



matter ( |Navarro, Frenk & White) [1997] ) and account for the effects of adiabatic 
contraction using a simple approach pioneered by Blumenthal et al. (1986), then 
one finds that for any WIMP mass less than 10 TeV, the effects of dark matter 
annihilation heating become significant during the collapse of the gas. | Spolyar et 



al. ( 2008 ) identify the point at which this occurs by comparing the heating rate due 



to dark matter annihilation with the H2 cooling rate. To determine a value for the 
latter, they make use of the simulation results of Yosh ida et al.| ([2006 ) and jGao et al.| 
(2007 ) and measure how the H2 cooling rate of the gas in the central collapsing core 
evolves as the collapse proceeds. They show that for a 100 GeV WIMP, heating 
dominates at gas densities n > 10 13 cm~ 3 . Finally, they argue that once dark matter 
annihilation heating dominates over H2 cooling, the gravitational collapse of the 
gas will come to a halt, and hence the gas will never reach protostellar densities. 
Instead, it will remain quasi-statically supported at a density of roughly 10 13 cm -3 
(for a 100 GeV WIMP), with a corresponding size scale of 17 AU, for as long as 
the dark matter annihilation rate remains large compared to the H2 cooling rate. As 
the time required to consume all of the dark matter within a radius of 17 AU may be 
hundreds of millions of years, the resulting quasi-static gas distribution - dubbed a 
"dark star" by Spolyar et al. ( |2008 1 - could potentially survive for a very long time. 
One criticism of the original |Spolyar et al. ( 2008[ ) model is its reliance on the 



Blumenthal et al. ( 1986) prescription for describing the effects of the adiabatic con- 



traction of the dark matter. This prescription assumes that all of the dark matter 
particles move on circular orbits, which is unlikely to be the case in a realistic dark 
matter minihalo, and concerns have been expressed that it may yield values for 
the dark matter density after adiabatic contraction that are significantly higher than 
the true values (see e.g. Gnedin et al. 2004) . For this reason, Freese et al. (2009 ) 
re-examined this issue using an alternative method for estimating the effects of adia- 
batic contraction, based on |Young] ( |1980| ). This alternative prescription does account 
for particles moving on radial orbits, and |Freese et al.| ( [2009| ) show that it predicts 
dark matter densities that are indeed systematically smaller than those predicted by 



the Blumenthal et al. ( 1986) prescription, but only by a factor of two. Freese et al. 



(2009 ) therefore conclude that although usi ng the |Young| (|1980|) pre scription for 
adiabatic contraction in place of the simpler Blumenthal et al. ( 1986) prescription 
will lead to some minor quantitative changes in the predicted outcome, the main 
qualitative results of the Spolyar et al. ( 2008| ) study are insensitive to this change, 
and one would still expect a "dark star" to form. 

Another potential problem with the dark star hypothesis is the fact that it is not 
at all clear that the collapse of the gas will stop once the dark matter heating rate 
exceeds the H2 cooling rate. For one thing, the values for the H2 cooling rate used 
by Spolyar et al. (2008 ) do not account for the effects of the dark matter annihilation 
heating. If this leads to an increase in temperature, then this will also increase the 
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H2 cooling rate, allowing more of the energy produced by dark matter annihilation 
to be radiated away. It is therefore unlikely that the point in the collapse at which 
the dark matter annihilation heating rate exceeds the Spolyar et al. estimate for the 
H2 cooling rate is marked by any sharp jump in the temperature. Instead, we would 
expect to find a more gradual temperature increase, at least up until the point at 
which collisional dissociation of the H2 starts to occur. 

Once H2 begins to dissociate, this provides another outlet for the energy gen- 
erated by dark matter annihilation. Spolyar et al. estimate that for a 100 GeV 
WIMP, the power generated by dark matter annihilation within the central core is 
Ldm ~ 140 L , and the core mass is roughly 0.6 M e . The total energy stored within 
the core in the form of the binding energy of the H2 molecules is roughly 

Eh, bind = 4.48eV x 0.76 x °- 6M ® ~ 2.6 x 10 45 erg, (56) 

and the time required for dark matter annihilation to produce this much energy is 

fdi s = %^ d ^200yr. (57) 

^dm 

For comparison, the free-fall time at this point in the collapse is roughly 15 years. 
H2 dissociation will therefore allow the collapse of the gas to continue until either 
the dark matter heating rate becomes large enough to destroy the H2 in the core in 
much less than a dynamical time, or the compressional heating produced during the 
collapse becomes capable of doing the same job. In either case, it is likely that much 
higher core densities can be reached than was assumed in the Spolyar et al. study. 
A first attempt to hydrodynamically model the formation of a "dark star" while 



correctly accounting for these thermodynamical effects was made by Ripamonti et 



[aTT|(|2010p. They used the ID, spherically symmetric hydrodynamical code described 



in Ripamonti et al. ( 2002 1 to model the collapse of the gas up to densities of order 
10 ls cm -3 for a range of different WIMP masses between 1 GeV and 1 TeV. Adi- 
abatic contraction of the dark matter was modelled using the algorithm described 
in Gnedin et al.| ( 2004) , and the effects of the dark matter annihilation heating and 



ionization were self-consistently accounted for in the chemical and thermal model. 
Ripamonti et al. (2010[) show that even in the most extreme case that they study, 



the heating produced by the dark matter appears unable to halt the collapse for an 
extended period. After the dark matter heating rate exceeds the H2 cooling rate, dis- 
sociation of H2 in the core accounts for most of the "excess" energy not radiated 
away by the gas, allowing the collapse to continue. Once the H2 in the core is ex- 
hausted, the temperature rises steeply, very briefly halting the collapse. However, the 
temperature quickly becomes large enough to allow other cooling mechanisms (e.g. 
H~ bound-free transitions or Lyman-a emission from atomic hydrogen) to operate, 
allowing the collapse to restart. Ripam onti et al.| ( |2010| l do not find any evidence for 
the formation of a hydrostatically supported "dark star" up to the highest densities 
that they study. Confirmation of this result in a 3D treatment of the collapse would 
be very useful. 
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2.3 The role of magnetic fields 
2.3.1 Initial strength 

The majority of the work that has been done on modelling the formation of the first 
stars assumes that magnetic fields play no role in the process, either because no mag- 
netic field exists at that epoch, or because the strength of any field that does exist is 
too small to be significant. A number of mechanisms have been suggested that may 
generate magnetic seed fields during the inflationary epoch, the electroweak phase 
transition or the QCD phase transition (see Kandus, Kunze & Tsagas 2011 for 



a recent comprehensive review). Observational constraints (e.g. Barrow, Ferreira & 



Silk 1997 Schleicher, Banerjee & Klessen 2008) limit the strength of the magnetic 



field at the epoch of first star formation to no more than about 1 nG (in comoving 
units), but it is quite possible that any primordial seed field resulting from one of 
these processes will actually have a much smaller strength. 

An alternative source for magnetic fields within the first generation of star- 
forming minihalos is the so-called Biermann battery effect ( |Biermanri 19501. In 



a partially ionized gas in which the gradient of electron density does not perfectly 
align with the gradient of electron pressure, as can happen if there is a tempera- 
ture gradient that does not align with the pressure gradient, the magnetic induction 
equation takes the form 

= Vx(vxB) + — ^ (58) 

at n%e 

where B is the magnetic field, v is the velocity, n e is the electron density, p e is the 
electron pressure, and e is the charge on an electron. In the limit that B 0, the first 
term on the right-hand side of this equation also becomes zero, but the battery term 
does not. It can therefore act as the source of a magnetic field in a gas that is initially 
unmagnetized. An early investigation into the effectiveness of the Biermann battery 
during galaxy formation was made by Kulsrud et al. ( 1997| l, who considered the 



formation of massive galaxies and showed that the Biermann battery coul d generate 



Xu et al. 



a field of strength B ~ 10~ 21 G during their assembly. More recently, 
(2008) have simulated the action of the Biermann battery during the formation of 
one of the first star-forming minihalos, finding that it is able to generate initial field 
strengths of the order of 10~ 18 G during this process. 



2.3.2 Amplification 

The seed fields generated by the Biermann battery, or by other processes acting 
in the very early Universe can be significantly amplified by flux-freezing during 
the gravitational collapse of the gas. If the diffusive timescale associated with am- 
bipolar diffusion or Ohmic diffusion is long compared to the gravitational collapse 
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timescale, then the magnetic field will be "frozen" to the gas, and will be carried 
along with it when the gas collapses. 

In the optimal case of spherical collapse, perfect flux freezing implies that the 
field strength evolves with density as B °c p 2 / 3 , and hence the magnetic pressure 
Pmag — j%% evolves as p mag <* p 4 / 3 . In comparison, the thermal pressure p t herm 
evolves as ptherm x pT, and so if the temperature does not vary much during the col- 
lapse, the plasma /3 parameter, j3 = ptherm/Pmag, evolves as /3 p~ 1,/3 . Therefore, 
if the gas is initially dominated by thermal pressure rather than magnetic pressure, 
it will remain so during much of the collapse, as a large change in the density is 
necessary to significantly alter p\ In the case examined by Xu et al. (2008 1, the very 
small initial magnetic field strength means that j3 is initially very large, and remains 
so throughout the collapse, implying that the magnetic field never becomes dynam- 
ically significant. Moreover, even if we take an initial comoving field strength of 
1 nG, comparable to the observational upper limit, at the mean halo density, /3 ~ 10 4 
(assuming a halo formation redshift z — 20 and a virial temperature of 1000 K), and 
does not become of order unity until very late in the collapse. Furthermore, if the 
collapse of the gas is not spherical, whether because of the effects of gravitational 
forces, angular momentum, or the influence of the magnetic field itself, the amplifi- 
cation due to flux freezing and collapse will be less than in the spherical case (see 
e.g. |Machida et al. 2006 who find a somewhat shallower relationship in some of 
their models). 

Therefore, for magnetic fields to play an important role in Pop. Ill star forma- 
tion, they must either start with a field strength very close to the observational up- 
per limit, or we must invoke an amplification process that is much more effective 
than the amplification that occurs due to flux freezing and gravitational collapse. 
One obvious possibility is amplification via some kind of dynamo process, which 
could bring about exponential amplification of an initially small seed field. Of par- 
ticular interest is the small-scale turbulent dynamo ( |Kraichnan & Nagarajan| |l967; 
Kazantsev)|196 8 ; Kulsrud & Anderson 1992| l. This produces a magnetic field that 



has no mean flux on the largest scales but that can have substantial mean flux within 
smaller-scale subregions. The growth rate of the magnetic field due to the turbulent 
dynamo is closely related to the rate of turnover of the smallest eddies. If the mag- 
netic field is sufficiently small that it does not significantly affect the velocity field 
of the gas (the kinematic approximation), and if we assume that we are dealing with 
Kolmogorov turbulence, then Kulsrud & Zweibel (2008 ) show that the magnetic en- 
ergy density grows exponentially, and that after a single gravitational free-fall time 
it is amplified by a factor exp(Re'/ 2 ), where Re is the Reynolds number of the flow. 
If we assume that the driving scale of the turbulence is comparable to the size of the 
mini halo, and that the turbulent vel ocity is of the same order as the sound speed (see 
e.g. Abel, Bryan & Norman 2002 1 >, then Re ~ 10 4 -10 5 , implying that the magnetic 
field is amplified by an enormous factor during the collapse. In practice, the field 
will not be amplified by as much as this analysis suggests, as the kinematic approx- 
imation will break down once the magnetic energy density becomes comparable to 
the kinetic energy density on the scale of the smallest eddies. Nevertheless, this sim- 
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pie treatment implies that the turbulent dynamo can amplify the magnetic field to a 
strength at which it becomes dynamically important. 

Although the importance of dynamo processes during the formation of the first 



galaxies has been understood for a number of years (see e.g. Pudritz & Silk 1989 



Becketal. 1994 Kulsrudetal. 1997 ), they have attracted surprisingly little atten- 
tion in studies of primordial star formation. Over the past couple of years, however, 
this has begun to change, with several recent studies focussing on the growth of 



magnetic fields during the formation of the first stars. The first of these was Schle- 
her et al. (2010), who studied the effectiveness of the turbulent dynamo during 
gravitational collapse using a simple one-zone Lagrangian model for the collapsing 
gas. Their model assumes that turbulence is generated by gravitational collapse on 
a scale of the order of the Jeans length, and that on smaller scales , the tu rbulent ve- 
locity scales with the length-scale I as v tur b « /^. |Schleicher et aL (2010 1 study both 
Kolmogorov turbulence, with J3 = 1/3 and Burgers turbulence, with j3 = 1/2, and 
show that in both cases, amplification of a weak initial seed field occurs rapidly, and 
that the field reaches saturation on all but the largest scales at an early point during 
the collapse. Because [Schleicher et al.| ( |2010| ) did not solve directly for the fluid ve- 
locities, they were unable to model the approach to saturation directly. Instead, they 
simply followed Subramanian ( 1998| ) and assumed that the strength of the saturated 
field satisfies 



Rm cr 



(59) 



where Rm cr ~ 60 is a critical value of the magnetic Reynolds number, Rm = v tur b/ / r) 
(where t] is the resistivity), that must be exceeded in order for exponential growth 
of the field to occur ( Subramanian) 1 1 998) . 

The main weakness of the Schleicher et al. ( 2010[ ) study lies in the assumptions 
that it was forced to make about the nature of the turbulent velocity field. Therefore, 



in a follow-up study, Sur et al. (2010) used high-resolution adaptive mesh refine- 
ment simulations to directly follow the coupled evolution of the velocity field and 



the magnetic field within a 3D collapse model. For their initial conditions, Sur et al. 
(|2010|) took a super-critical Bonnor-Ebert sphere (Bonnor 1956||Ebert||1955[) with 



a core density n c 



10 4 cm~ 3 and a temperature T = 300 K. They included initial 
solid-body rotation, with a rotational energy that was 4% of the total gravitational 
energy, and a turbulent velocity component with an RMS velocity equal to the sound 
speed and with an energy spectrum £ (A:) °< k~ 2 . A weak magnetic field was also in- 
cluded, with an RMS field strength B rms = 1 nG, and with the same energy spectrum 
as the turbulence. For reasons of computational efficiency, Sur et al. ( 2010| ) did not 
follow the thermal and chemical evolution of the gas directly. Instead, they adopted 
a simple barotropic equation of state, P«< p 1 A , inspired by the results of previous 
hydrodynamical models (e.g. |Abel, Bryan & Nor man 2002 ). In view of the sensi- 
tivity of the turbulent dynamo to the Reynolds number, and the fact that numerical 
dissipation on the grid scale limits the size of Re in any 3D numerical simulation 
to be substantially less than the true physical value, there is good reason to expect 
that the dynamo amplification rate will be sensitive to the numerical resolution of 
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the simulation. |Sur et aT7| ( |201Q) > therefore focussed on the effects of resolution, per- 
forming five different simulations with the same initial conditions, but with different 
grid refinement criteria. Starting with a model in which the refinement criterion en- 
sures that the Jeans length is always resolved by 8 grid zones, they looked at the 
effects of increasing this number to 16, 32, 64 and 128 grid zones. 

|Sur et al.] (|2010) showed that in the 8 and 16 cell runs, the magnetic field strength 
increases with density at a slower rate than the B« p 2 / 3 that we would expect sim- 
ply from flux freezing and roughly spherical collapse, indicating that in these runs, 
the turbulent dynamo does not operate. Starting with the 32 cell run, however, they 
found evidence for an increase in B with density that is larger than can be explained 
simply by compression, which they ascribe to the effects of the turbulent dynamo. 
They showed that as the number of grid zones used to resolve the Jeans length is 
increased, the rate at which the field grows also increases, and there is no sign of 
convergence at even their highest resolution. This resolution dependence explains 
why the earlier study of Xu et al. ( 2008| l found no evidence for dynamo amplifica- 
tion, as their study used only 16 grid zones per Jeans length. 

More recently, Fed errath et al.| ( [2~01 \\ have re-examined this issue of resolution 
dependence, and have shown that when the number of grid zones per Jeans length 
is small, the amount of turbulent energy on small scales is also significantly under- 
estimated. The reason for this is the same as the reason for the non-operation of 



the turbulent dynamo: the effective Reynolds number is too small. Federrath et al. 



( |201 l| l show that in gravitationally collapsing regions that undergo adaptive mesh 
refinement, the effective Reynolds number scales with the number of grid zones 
per Jeans length as Re e ff = (N/2) 4 / 3 . Furthermore, an effective Reynolds number 
Re e ff ~ 40 is required in order for the turbulent dynamo to operate, implying that 
one needs ^ 30 or more zones per Jeans length in order to begin resolving it, in 
agreement with the findings of |Sur et~aL ( 2010| >. It should also be noted that the 
operation of the turbulent dynamo in simulations of turbulence without self-gravity 



requires a similar minimum value for the Reynolds number (Haugen, Brandenburg 
|& Dobler||2004l >. 

Together, these studies support the view that amplification of a weak initial mag- 
netic field by the turbulent dynamo may indeed have occurred within the first star- 
forming minihalos. However, a number of important issues remain to be addressed. 
First, the three-dimensional studies carried out so far all adopt a simple barotropic 
equation of state, rather than solving self-consistently for the thermal evolution of 
the gas. This is a useful simplifying assumption, but may lead to incorrect dynami- 
cal behaviour, as one misses any effects due to thermal instabilities, or the thermal 
inertia of the gas (i.e. the fact that the cooling time is typically comparable in size 
to the dynamical time). Work is currently in progress to re-run some of these mod- 
els with a more realistic treatment of the thermodynamics and chemistry in order 
to explore the effect that this has on the degree of amplification (T. Peters, private 
communication). Second, it will clearly be important to perform similar studies us- 
ing more realistic initial conditions for the gas. Of particular concern is whether the 
turbulence that is generated during the gravitational collapse of gas within a primor- 
dial minihalo is similar in nature to that studied in these more idealized calculations, 
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and if not, what influence this has on the amplification of the field. Finally, and most 
importantly, there is the issue of the level at which the field saturates. Exponential 
amplification of the field by the small-scale dynamo will occur only while the kine- 
matic approximation holds, i.e. while the energy stored in the magnetic field is much 
smaller than the energy stored in the small-scale turbulent motions. Once the field 
becomes large, the Lorentz force that it exerts on the gas will act to resist further 
folding and amplification of the field. In addition, the dissipation of magnetic en- 
ergy by Ohmic diffusion and ambipolar diffusion will grow increasingly important. 
However, it remains unclear which of these effects will be the most important for 
limiting the growth of the magnetic field in dense primordial gas. 



2.3.3 Consequences 

If a strong magnetic field can be generated by dynamo amplification, then it will 
affect both the thermal and the dynamical evolution of the gas. The possible dy- 
namical effects of a strong magnetic field have been investigated by Ma chida et al.| 



(2006 2008). In a preliminary study, Machida et al. ( 2006) 1 used nested-grid sim 



ulations to investigate the influence of a magnetic field on the collapse of a small, 
slowly-rotating primordial gas cloud. For their initial conditions, they used a super- 
critical Bonnor-Ebert sphere with mass 5.1 x 1O 4 M0, radius 6.6 pc, central density 
n c = 10 3 cm~ 3 and an initial temperature of 250 K. They assumed that this cloud 
was in solid body rotation with angular velocity Qq and that it was threaded by 
a uniform magnetic field oriented parallel to the rotation axis, with an initial field 
strength Bo. They performed simulations with several different values of £2q and 
Z?o, with the former ranging from 10~ 17 s _1 to 3.3 x 10~ 16 s -1 , and the latter from 
10~ 9 G to 10~ 6 G. To treat the thermal evolution of the gas, they used a barotropic 
equation of state, based on the one-zone results of |Omukai et al.| ( |2005) . 

Mach ida et al.|p006| used this numerical setup to follow the collapse of the gas 



down to scales of the order of the protostellar radius. They showed that the mag 
netic field was significantly amplified by compression and flux freezing during the 
collapse, reaching strengths of order 6 x 10 5 -6 x 10 6 G on the scale of the proto- 
star. A very compact disk with a radius of few R E formed around the protostar, and 
in models with initial field strength Bq > 10~ 9 G, the magnetic field became strong 
enough to drive a hydromagnetic disk wind that ejected roughly 10% of the infalling 
gas. Numerical limitations (discussed in Section |3]below) prevented Machida et al. 



( 2006 ) from following the evolution of the system for longer than a few days after 



the formation of the protostar, and so it remains unclear whether an outflow would 
eventually be generated in the 10~ 9 G case, and whether the outflows continue to be 
driven as the protostar and disk both grow to much larger masses. 

In a follow-up study, |Machida et al.| ( |2008] > used a similar numerical setup, but 
examined a much broader range of values for j3o(= £rot/|£grav|), the ratio of the 
initial rotational energy to the initial gravitational energy, and 7o(= E mag /\E gav \), 
the ratio of the initial magnetic energy to the initial gravitational energy. They found 
that the outcomes of the simulations could be classified into two main groupings. 
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Clouds with j3o > 70, i.e. ones which were rotationally dominated, formed a promi- 
nent disk during the collapse that then fragmented into a binary or higher order 
multiple system. In these simulations, no jets were seen (with the exception of a 
couple of model in which j3rj ~ 7o). On the other hand, when /3o < 7o, i.e. when the 
cloud was magnetically dominated, the disk that formed was much less prominent 
and did not fragment, but instead an outflow was driven that removed of order 10% 
of the gas that reached the disk, as in the Machida et al. (2006 ) study. 

These results support the idea that outflows will be a natural consequence of the 
generation of strong magnetic fields during Population III star formation. However, 
it is important to note that the Machida et al. ( 2006 , 2008 ) simulations only model 
the very earliest stages of outflow driving, on a timescale t -C 1 yr. The evolution of 
outflows on much longer timescales, and their influence on the infalling gas have not 
yet been studied in detail, and it is unclear to what extent one can safely extrapolate 
from the very limited period that has been studied. 

A strong magnetic field could also have a direct impact on the thermal evolution 
of the gas, through the heating arising from ambipolar diffusion. The effects of this 
process in gravitationally collapsing gas within the first star-forming minihalos have 
been investigated by Schleicher et al. ( 2009} using a simple one-zone treatment of 
the gas. They assume that in the absence of ambipolar diffusion, the magnetic field 



strength would evolve as5<x p a , where a = 0.57(Mj/Mj ma < 



\0.0116 



and Mj, ma{ , is 



the magnetic Jeans mass (i.e. the minimum mass that a perturbation must have in 
order to be unstable to its own self-gravity when support against collapse is pro- 
vided by a magnetic field, rather than by thermal pressure). This expression for a 
is an empirical fitting formula derived by Schleicher et al. ( 2009 1 from the results 
of Machi da et al. ( 2006| l. In their treatment of the evolution of B within their one- 
zone models, Schleicher et al. ( 2009| l also account for the loss of magnetic energy 
through ambipolar diffusion. 

Another important simplification made in the Schleic her et al.| ( |2009| l model is 



the replacement of the full expression for the ambipolar diffusion heating rate (Pinto 

[eTaTj [20081 
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and tjad is the ambipolar diffusion resistivity, with the simpler ap- 
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where Lg is the coherence length of the magnetic field. 



(61) 



Schleicher et al. (2009 1 show that if the initial field strength is less than 0.01 nG 



(comoving), then ambipolar diffusion heating has almost no effect on the thermal 
evolution of the gas. For stronger fields, ambipolar diffusion heating leads to an in- 
crease in the gas temperature at densities between n ~ 10 4 cm~ 3 and n ~ 10 10 cm -3 , 
amounting to as much as a factor of three at n ~ 10 8 cm~ 3 . However, at higher 
densities, three-body H2 formation heating becomes a more important heat source 
than ambipolar diffusion, meaning that the temperature evolution becomes largely 
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independent of the magnetic field strength once again. Schlei cher et al.| ( |2009| l did 
not examine initial field strengths larger than 1 comoving nG, as these are ruled 
out by observational constraints, but if one considers the effects of the turbulent dy- 
namo acting during the collapse, then it is possible that much larger fields could be 
generated on smaller scales, and it would be useful to revisit this issue and examine 
whether ambipolar diffusion heating from these smaller-scale fields can significantly 
affect the collapse of the gas. 

Finally, one important caveat to bear in mind regarding the Machida et al. ( 2006 , 
2008) and Schleicher et al. (2009) simulations is the fact that they adopt a correlated 



initial magnetic field, while the field generated by the turbulent dynamo will have 
little or no correlation on large scales (Maron, Cowley & Mc Williams , 2004). The 
extent to which the dynamical and thermal effects of this uncorrected field will be 
the same as those of a correlated field is unclear. Further investigation of this issue 
would be very valuable. 



3 Evolution after the formation of the first protostar 

As we saw in the last section, when it comes to the formation of the very first Pop- 
ulation III protostar, there is broad agreement on the details of the process, with 
different groups, who use different numerical approaches, finding results that are in 
good qualitative agreement with each other. Some quantitative disagreements still 



exist (see e.g. Turk et al. 2011 1, but it is unclear to what extent these reflect real 
differences between numerical approaches as opposed to natural variation in the de- 
tails of the collapse. The main uncertainties in this phase stem from uncertainties 
in the input physics, such as whether magnetic fields can become amplified to dy- 
namically significant levels during the collapse, or whether dark matter annihilation 
significantly affects the outcome. 

Once we move on to considering the evolution of the gas within star-forming 
minihalos after the formation of the first protostar, the situation becomes much less 
clear. The fundamental problem stems from the fact that although we can follow the 
gravitational collapse of the primordial gas down to scales as small as the protostel- 
lar radius (see e.g. Yoshida, Omukai & Hernquist 2008), the numerical timestep in 



an explicit hydrodynamical code becomes extremely short during this process. This 
is a consequence of the Courant condition, which states that for such a code to be 
numerically stable, the timestep must satisfy 

Ax 

At<— , (62) 
c s 

where Ax is the size of the smallest resolution element, and c s is the sound speed of 
the gas. 

The Courant condition implies that if we take a value of Ax small enough to ade- 
quately resolve the structure of the protostar and the gas immediately surrounding it 
(e.g. Ax = 1R©), then the required timestep will be extremely small: At < 7 x 10 4 s 
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for Ax = lR e and a sound speed of lOkms . This means that if we want to fol- 
low the later evolution of the protostar and the surrounding gas over a timescale of 
thousands of years in order to see how it grows in mass prior to reaching the main 
sequence, then we must use a very large number of timesteps: our simple estimate 
above yields a number of the order of a million. In practice, the computational ex- 
pense of doing this within a three-dimensional hydrodynamical code is prohibitively 
large, meaning that it has so far proved impossible to study the evolution of the gas 
in this fashion. 

Efforts to surmount this difficulty typically follow one of two approaches. One 
approach is simply to halt the numerical simulation at the point at which the Courant 
timestep becomes prohibitively small, and to model the later evolution of the pro- 
tostar using a semi-analytical, or one-dimensional, fully numerical treatment. To do 
this, it is necessary to make some assumption about the behaviour of the gas sur- 
rounding the protostar. In general, models of this type assume that the gas does not 
fragment and form additional protostars, but instead is simply accreted by the exist- 
ing protostar, either directly or via a protostellar accretion disk. The results obtained 
using this approach - what we afterwards refer to as the "smooth accretion model" 
- are discussed in Section IXTl below. 

The other approach that can be used to study the further evolution of the gas sur- 
rounding the protostar makes use of a technique developed for studies of contem- 
porary star formation, which face a similar problem on protostellar scales. Gravita- 
tionally bound regions of gas that become smaller than some pre-selected size scale 
are replaced by what are usually termed sink particles (see e.g. |Bate, Bonnell & 
Price|[T995] l. These particles can accrete gas from their surroundings and continue 
to interact gravitationally with the surrounding gas, but allow one to neglect the very 
small-scale hydrodynamical flows that would otherwise force one to take very small 
numerical timesteps owing to the Courant condition. The great advantage of the sink 
particle technique is that one need make no assumption about the dynamical evolu- 
tion of the gas surrounding the protostar on scales much larger than the effective size 
of the sink particle (the so-called accretion radius, discussed in more detail below), 
as one can simply continue to model this using the same numerical techniques as 
were used to model the initial gravitational collapse. The main disadvantage of the 
technique is that, strictly speaking, it represents an ad hoc modification of the fluid 
equations, with consequences that may not be entirely straightforward to predict. 
The modification to the solution caused by replacing dense gas with sink particles is 
unlikely to significantly affect the evolution of the gas on scales that are much larger 
than the accretion radius, but will clearly have an effect on the flow on scales close 
to the accretion radius. In addition, the common strategy of treating sink particles as 
point masses may not be appropriate when dealing with close encounters between 
sinks, as one misses the tidal forces acting between the gas clumps represented by 
the sinks. 

Although sink particles have been used in studies of Population III star formation 
for over a decade, simulations using the correct initial conditions, and with sufficient 
spatial resolution and mass resolution to capture the details of the gas flow on scales 
close to those of individual protostars have only recently become possible. These 
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simulations show that, contrary to the assumption made in the smooth accretion 
model, the gas generally fragments, rather than simply accreting onto a single, cen- 
tral protostar. The results obtained from studies using sink particles - afterwards 
referred to as the "fragmentation model" for Population III star formation - are dis- 
cussed in Section [3~2l below. 



3.1 The smooth accretion model 



3.1.1 Determining the accretion rate 



As we have already discussed above, at the point at which the protostar forms, its 
mass is very small (M ~ 0.01 M ; see Yoshida, Omukai & Hernquist|20 08 ), but it 



is surrounded by an infalling envelope of gas containing tens or hundreds of solar 
masses. If we assume that the gas in this infalling envelope does not undergo gravi- 
tational fragmentation, then it has only two possible fates - it must either be accreted 
by the central protostar (or protostellar binary; see e.g. Turk, Abel & O'Shea 2009) 



or it must be prevented from accreting, and possibly expelled from the immediate 
vicinity of the protostar, by some form of protostellar feedback. This means that the 
mass of the protostar at the point at which it forms has very little to do with its final 
mass. To determine the size of the latter, we must understand the rate at which gas is 
accreted by the protostar, and how this process is affected by protostellar feedback. 

Since protostellar feedback involves a number of different processes, many of 
which are complicated to model, it is easiest to start by considering models in which 
feedback effects are not included. As feedback acts to reduce the accretion rate, 
models of this type allow us to place an upper limit on the final mass of the Pop. Ill 
star. 

A useful starting point is a simple dimensional analysis. Suppose that the proto- 
star is embedded in a gravitationally unstable cloud of mass M and mean density 
(p), and that the protostellar mass -C M, so that its gravity is negligible in com- 
parison to the self-gravity of the cloud. The timescale on which the gas cloud will 
undergo gravitational collapse and be accreted by the protostar is simply the free- 
fall collapse time, fff = y / 3n/32G(p). Therefore, the time-averaged accretion rate 
will be given approximately by 

M / 

M M ~ — ~My/G{pj. (63) 
fff 

If the gas cloud were highly gravitationally unstable, then it would fragment rather 
than accreting onto a single object, so let us assume that it is only marginally un- 
stable, i.e. that M ~ Mj. In that case, since Mj ~ c^G~ 3 / 2 p~ 1//2 , we can write our 
estimate of the time-averaged accretion rate as 



(64) 
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(65) 



We therefore find that the characteristic accretion rate scales as the cube of the sound 
speed. Moreover, since c s °^ T 1 / 2 , this implies that the accretion rate scales with 
temperature as M °^ 

This is an important result, because as we have already seen, the characteristic 
temperature of the dense, star-forming gas in a primordial minihalo is of the order 
of 1000 K, far larger than the 10 K temperatures found within prestellar cores in 
local regions of star formation (see e.g. Bergin & Tafalla, 2007). Our simple scaling 
argument therefore tells us that we will be dealing with far higher accretion rates in 
the Population III case than we are used to from studies of local star formation. 

If we want to improve on this simple scaling argument and derive a more accu- 
rate figure for the accretion rate, then there are three main ways in which we can go 
about it. One possible approach is to construct a simplified model for the collaps- 
ing protostellar core from which an approximation to the true accretion rate can be 
derived analytically (or with only minor use of numerical calculations). For exam- 
ple, if we assume that the protostellar core is isothermal and spherically symmetric, 
then there is a whole family of similarity solutions that could potentially be used to 
describe the collapse (Hunter 1977 Whitworth & Summers| fl985), including the 
familiar Larson-Penston solution~ ^Larson| 1 1*969 ; Penston, 1969), or the Shu solution 
( |Shul[T977l >. 



An example of this approach is given in Omukai & Nishi (1998). These au- 
thors used a spherically symmetric Lagrangian hydrodynamical code to simulate 
the formation of a Population III protostar, and found that prior to core formation, 
the gravitational collapse of the gas could be well described with a Larson-Penston 
similarity solution, with an entropy parameter K = p/p? — 4.2 x 10 1 1 (in cgs units) 
and an effective adiabatic index y e ff = 1.09. Omukai & Nishij ( |T998 1 were unable 
to continue their numerical study past the point at which the protostar formed, for 
the reasons addressed above, but assumed that the same similarity solution would 
continue to apply. By making this assumption, they were therefore able to derive the 
following accretion rate for the protostar 



M = 8.3 x 10" 



t 

1 yr 



-0.27 



M yr" 



(66) 



M = 6.0x 10" 



-0.343 



In a similar study, using a more sophisticated treatment of the microphysics of the 
collapsing gas, iKip amonti et al.| ( |2002| ) also found that the initial flow was well 
described as a Larson-Penston similarity solution, but derived a different accretion 
rate 

,(± 

Another example of this approach comes from Tan & McKee ( 2004 ). They model 
the accretion flow onto a Pop. Ill protostar as a spherical, isentropic polytrope, and 
derive an accretion rate that is a function of three parameters: the entropy parameter 



M yr" 



(67) 
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K, the polytropic index y p (which, for an isentropic flow, is equal to the adiabatic 
index y), and 0*, a numeric parameter of order unity, which is related to the initial 
conditions of the flow. Tan & Mc Kee (2004) use the results of Omukai & Nishi 



( 1998 1 and Ripamonti et al. 1(12002 1 to argue that y p = 1.1, and use the 3D simulation 



2002 ) to set the other two parameters to 0* = 1 .43 



results of Abel, Bryan & Norman 
and K — 1.88 x \0 vl K' (in cgs units), where 

and where the effective temperature T e ff = P e s/(nk) accounts for the contribution 
made to the pressure by small-scale, subsonic turbulence in addition to the standard 
thermal pressure. Based on this, they then derive the following rate for the accretion 
of gas onto the protostar and its associated accretion disk 

/ , \ -°- 30 



M = 7.0 x W^K ' I — J Moyr" 1 . (69) 

This can be directly compared to the other determinations of M if we assume that 
all of the gas reaching the accretion disk is eventually accreted by the star, which is 
a reasonable assumption for models that do not include the effects of gravitational 
fragmentation or protostellar feedback. 

Instead of using simulation results to select a particular collapse model (e.g. 
Larson-Penston collapse) and then calculating M from the model, the second main 
approach used to determine M attempts to infer it from the state of the gas in the sim- 
ulation at the point at which the protostar forms, using the information that the sim- 
ulation provides on the density and velocity distributions of the gas. This approach 
was pioneered by |Abel, Bryan & N orman (2002), who considered two simple mod- 
els for the time taken for a given fluid element to accrete onto the central protostar. 
In the first of these models, they assumed that the time taken for the gas within a 
spherically-averaged shell of radius r to accrete onto the protostar was given by the 
ratio between the mass enclosed within the shell, M(r), and the rate at which gas 
was flowing inward at that radius, i.e. 

• - MW (70) 



Anr 2 p{r)\v r {r)[ 



where p(r) and v r (r) are the spherically-averaged density and radial velocity in the 
shell. In the second model for f acc , they used an even simpler approximation, setting 
4cc to the time that it would take for the gas to reach the protostar if it merely 
maintained its current radial velocity, i.e. 



face--- (71) 

Vr 



Abel, Bryan & Norman ( 2002) show that other than at the very earliest times, these 



two approaches yield very similar values for f acc , and hence very similar values 
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for the accretion rate. This strategy has subsequently been used by many other au- 
thors to derive predicted protostellar accretion rates from their simulations (see e.g. 



Yosh ida et al^[2006l[O 7 Shea & Norman] [2007) |McGreer & Bryan1|2TX)8l|TurkeUL 
201 1] >. Of particular note is the study by O'Shea & Norman| ( |2007| l,~who perform 
multiple simulations of Population III star formation using different random real- 
izations of the cosmological density field. They find that minihalos assembling at 
higher redshifts form more H2 than those assembling at lower redshifts, owing to 
the higher mean density of the virialized gas in the high redshift minihalos. They 
show that in their simulations, this leads to the gas at densities n > 10 4 cm~ 3 having 
significant differences in its mean temperature in the different halos. In the most 
H2-rich minihalos, the dense gas can be as cold as 200 K, while in the minihalos 
with the least H2, it can be as high as 1000 K. As a result, the predicted accretion 
rates for the different minihalos span more than an order of magnitude, thanks to the 
strong scaling of M with temperature. Unfortunately, it is necessary to treat these re- 



sults with a degree of caution, as the O'Shea & Norman (2007) simulations did not 
include the effect of three-body H2 formation heating, which is known to have a 
significant influence on the temperature of the dense gas. It is unclear whether sim- 
ulations that include this effect would produce dense gas with such a wide range 
of temperatures and accretion rates, although a study that is currently being carried 
out by Turk and collaborators should address this issue in the near future (M. Turk, 
private communication). 

The third main approach used to determine the accretion rate involves measuring 
it directly in a simulation of the later evolution of the gas around the protostar. If 
we replace the protostar with a sink particle, then we can measure M simply by 
measuring the rate at which the sink particle mass increases. This approach was 
first used by Bromm & Loeb ( 2004|, in a study of Population III star formation in 
which a sink particle was created once the gas density exceeded a threshold value 



«th = 10 12 cm 3 (we will have more to say about this study below). Bromm & Loeb 



( 2004 ) showed that the rate at which gas was accreted by the sink particle could be 
approximated as a broken power-law 



3.6 10 2 f ^ 



-0.25 



M yr 



-1 



M = 



6.3 x 10- 1 MQyr- 1 



t < 10 3 yr 



t > 10 3 yr 



(72) 



for times t < 10 4 yr. Bromm & Loeb halted their simulation at t ~ 10 4 yr and hence 
could not directly measure the evolution of M at later times, although they did con- 
sider what the final mass of the protostar would be if one simply extrapolated Equa- 
tion r 



72 over the three million year lifetime of a massive star. 



Accretion rates have also been measured using the sink particle technique in the 
group of simulations carried out by |Clark et al.| l |2011a|b) , |Greif et aL] ( |201 1 a] > and 
Smith et al. ( 201 1 J > that find evidence for fragmentation of the gas (see Section 3.2 
below). The accretion rates onto the individual sinks show a considerable degree of 
variability in these calculations, but the total accretion rate, i.e. the rate of change of 
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the sum of all of the sink particle masses, evolves more smoothly with time, and is 
of a similar order of magnitude to the other estimates plotted above. 




Time (yr) 



Fig. 3 Three different estimates for the accretion rate onto a Pop. Ill protostar, taken from Tan & 
McKee (2004; solid line), Turk et al. (2011; dashed line) and Smith et al. (2011; dotted line), as 
described in the text. Results from the Smith et al. 1 201 1 1 simulation are only plotted for the period 
covered by the simulation, i.e. t < 2 x 10 4 yr. 



In Figure[3]we compare several of these different estimates for the accretion rate. 



We plot three examples, derived using different techniques: a rate based on the Tan 
& McKee] ( |2004| l formalism, computed assuming that K' = 1; a rate inferred from 



the results of one of the adaptive mesh refinement simulations presented in Turk 



et al. ( 201 l|l - specifically, the simulation that was run using the |Palla, Salpeter, & 



Stahler (1983) rate coefficient for three-body H2 formation; and a rate measured 



using sink particles, taken from Smith et al. ( 201 l| l. 

At very early times (t < lOOyr), the three different techniques yield rather differ- 
ent estimates for M, but this is primarily a consequence of the limited resolution of 
the numerical simulations. At later times, we see that both the Tan & McKee (2004 ) 
formalism and the |Turk et al. ( 2011| l simulation predict a similar form for the ac- 
cretion rate, but disagree by about a factor of two on the normalization, which may 
simply indicate that our adopted value of K' is slightly too large. We also see the 
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same general trend in the |Smith et al.| ( [201 1| > results, but in this case there is consid- 
erable and rapid variation in M with time. This is a result of the fragmentation of the 
gas in this simulation, which produces a set of sink particles that undergo chaotic 



N-body interactions (see Section 3.2 below). A similar effect is seen in simulations 



of protostellar accretion in present-day star-forming regions (see e.g. Stamatellos, 
IWhitworth & Hubber||20TT) . 

Regardless of whether M varies smoothly or erratically with time, one fact that 
is clear from Figure [3] is that the protostellar accretion rate remains very large for a 
considerable time. This implies that the total mass of gas that is converted to stars 



can become fairly large after a relatively short time. For example, if we take the Tan 



& McKee ( 2004 1 estimate with K' = 1 as a guide, then we find that the total mass 



, 0.70 



in stars increases with time as: 



M*=0.1I^-) Meyr 1 . (73) 

This means that after 5 x 10 4 yr (the Kelvin-Helmholtz relaxation time for a 1OOM 
star), we have ~ 195 M , while after 2 x 1 6 yr (the typical lifetime for an O star), 
we have ~ 2575 M . Therefore, if the gas does not fragment and protostellar 
feedback is ineffective, one is led to the prediction that the resulting Population III 
star will be extremely massive. In practice, the gas probably does fragment (see 



Section 3.2 below), and protostellar feedback cannot be completely ignored, but 
even so, we would expect to be able to form massive Population III stars relatively 
easily. 

Finally, it should be noted that so far we have considered accretion only in the 
standard ^-dominated case, i.e. in a minihalo with a minimum gas temperature of 
around 200 K. In minihalos that reach much lower temperatures through HD cool- 
ing, the predicted accretion rates are smaller, as one would expect from the simple 
dimensional analysis given at the start of this section. For example, if one uses the 



Tan & McKee (2004) formalism to estimate the accretion rate, then Equation 69 
still applies, but the value of K' is significantly smaller. Taking n = 10 6 cm~ 3 and 



r e ff = 150 K as plausible values to substitute into Equation 68 we find that K' ~ 0.3, 



and hence the predicted accretion rate is roughly a factor of six smaller than in the 



H2-dominated case. Values estimated from numerical simulations using the Abel 



Bryan & Norman| d2002 ) approach agree fairly well with this simple estimate (see 



e.g. Yoshida, Omukai & Hernquistl|2007|pcGreer & Bryan 2008) 



3.1.2 Protostellar structure and evolution 

Having established how quickly gas will be accreted by the protostar in the ab- 
sence of feedback, the obvious next step is to examine how this will be modified by 
protostellar feedback. Before doing this, however, we must first spend a little time 
discussing what is known about the internal structure of Population III protostars, 
and how this evolves with time. This is important if we want to understand how the 
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radius and luminosity of a given Pop. Ill protostar evolve, and these quantities are 
obviously of great importance when determining the influence of that protostar on 
the surrounding gas. 

The internal structure of a Pop. Ill protostar, and how this evolves as the protostar 



ages and accretes matter from its surroundings was first studied in detail by Startler, 
Palla & Salpeter ( 1986a b). They assume that the accretion process can be treated 
as a series of quasi-steady-state accretion flows onto a hydrostatic core, which is 
bounded by a strongly radiating accretion shock. Within the core, the standard stel- 
lar structure equations are solved. Outside of the core, the treatment depends on the 
optical depth of the gas. If the gas is optically thin to the radiation from the accre- 
tion shock, then the accretion flow is assumed to be in free-fall. Otherwise, a more 
detailed calculation is made that incorporates the effects of the radiation force on 
the infalling gas. The accretion shock itself is treated as a simple discontinuity. 



In their initial study, Stahler, Palla & Salpeter ( 1986a i began with a core mass of 
0.01 M and followed the growth of the protostar until its mass reached 10.5 M . 
They assumed a constant accretion rate M = 4.41 x 10~ 3 M Q yr _1 , and found that 
for this choice of accretion rate, the evolution of the protostar could be divided into 
three qualitatively distinct phases. 

In the first phase, which lasts until the protostellar mass M* =0.1 M , the pro- 
tostar relaxes from its initial entropy profile into one consistent with the selected 
Stahler, Palla & Salpeter| ([1986a) dub this a 'decay of transients' 



accretion rate. 



phase, and the fact that it quickly comes to an end shows that although the initial 
conditions used in the Stahler, Palla & Salpeter ( 1986a| l study are probably incorrect 



in detail, the flow soon loses all memory of them, and therefore any inaccuracy at 
this stage is unlikely to affect the later results. 

Once the initial transients have died away, the protostar enters the second phase of 
its evolution. During this phase, its central temperature remains low(r c -10 5 K), re- 
sulting in a high interior opacity and hence a low interior luminosity. Consequently, 
the evolution of the core during this phase is almost adiabatic; although the core 
continues to gradually contract, this contraction does not lead to any increase in the 
central entropy. Since the postshock entropy increases over time due to the increas- 
ing strength of the accretion shock (which is itself a natural result of the increasing 
protostellar mass), the core develops an off-centre distribution of entropy and tem- 
perature. 

The gas surrounding the accretion shock remains optically thick throughout this 
period. This is a direct result of the high accretion rate, which produces a highly 
luminous accretion shock. This produces sufficient radiation to partially ionize the 
preshock gas in the vicinity of the shock, creating a structure known as a radiative 
precursor. The H opacity of the dense, partially ionized gas in this radiative pre- 



cursor is more than sufficient to make it optically thick. Stahler, Palla & Salpeter 
show that the core radius during this period evolves as 



R * = 48A [yi) Uo» R - (74) 
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where Mq = 4.41 x 10 M yr , while the photospheric radius evolves as 

M t \ - 27 /M\ 0M 



^ = 66 ' 8 UU UoJ R0 ' (75) 

so Rp > throughout. The strong H opacity also keeps the photospheric tem- 
perature low (T p ~ 6000 K), which prevents the protostar from being able to ionize 
material outside of its photosphere. 

This near-adiabatic accretion phase comes to an end once the cooling time of the 
core, given approximately by the Kelvin-Helmholtz timescale 

*e-j£, (76) 

becomes comparable to the accretion timescale f acc = M*/M. This occurs for a core 
mass M ~ 1M B , and results in the core entering a phase of homologous collapse, 
while energy and entropy are transferred outwards in the form of a 'luminosity 
wave' . The radial position of the luminosity peak moves outwards towards the accre- 
tion shock, reaching it at about the time that the core mass has reached 8M . This 
results in a rapid swelling of the outermost layers, which weakens the accretion 
shock and leads to it becoming optically thin. Stahler, Palla & Salpeter terminate 
their simulation shortly afterwards, once the core mass has reached 10. 5M©. 

Stahler, Palla & Salpeter ( 1986b| l simulate the later stages of the evolution of a 



primordial protostar. Their initial protostellar core has a mass of 5M Q , and they 
evolve this core forward in time, assuming that no further accretion occurs (i.e. the 
protostellar mass remains fixed at 5M ). They find that deuterium burning within 
the protostar begins after only 6000 years, but that hydrogen ignition does not occur 
until f = 2 x 10 5 yr, and the protostar does not reach the zero-age main sequence 
(ZAMS) until t ~ 10 6 yr. 

An improved treatment of the later stages of the evolution of the protostar 



was made by |Omukai & Paha (2001 1. They used a very similar setup to that in 
Stahler, Palla & Salpeter (1986a), albeit with improved zero metallicity opacities, 
and adopted the same constant accretion rate, M = 4.41 x 10~ 3 M© yr~ 1 . However, 



unlike Stahler, Palla & Salpeter ( 1986a), they initialized their simulation at the point 



at which the core mass was M = 8 M©, but did not halt the simulation once the core 
had grown to 10.5 M©. Instead, they continued to follow the growth of the protostar 
until well after hydrogen ignition. They found that deuterium burning within the 
core began once the core mass was 12M (corresponding to a time t = lOOOyr after 
the beginning of the simulation, given the assumed accretion rate), and that it was 
complete by the time the mass had reached 30 M (corresponding to t = 5000 yr). 
Hydrogen ignition followed roughly 11000 years later, at t — 1.6 x 10 4 yr after the 
beginning of the simulation, at which time the mass of the protostar was 80 M©. 
At this point, the internal luminosity of the protostar is very close to the Eddington 
value, which leads to the outer layers of the protostar developing oscillatory be- 
haviour: the high luminosity leads to expansion, the expansion causes the accretion 
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luminosity to drop, the reduced luminosity can no longer maintain the expansion, 
leading to contraction of the core, and the contraction raises the accretion luminos- 
ity, allowing the whole cycle to begin again. Finally, once the core mass reaches 
300 M Q , at t ~ 6.6 x 10 4 yr, the contribution of nuclear burning to the protostellar 
luminosity becomes large enough to drive a final phase of expansion that is strong 
enough to terminate accretion onto the protostar. Omukai & Palla (2001 ) halt their 
simulation at this point. 



In a follow-up study using a similar spherically-symmetric setup, Omukai & 



Palla 



( 2003 i performed the same analysis for a range of different values of M, look- 
ing at models withM = (0.25,0.5, 1.0,2.0) x M M (where M fid was the rate adopted 
by |Stahler, Palla & Salpeter|1986a| and |Omukai & Palla|2001| >, as well as a model 
using the time-dependent accretion rate predicted byjAbel, Bryan & Norman (2 002[ ). 
The earliest stages of protostellar evolution are qualitatively the same in all of these 
models: we see again the same sequence of adiabatic growth, propagation of a lumi- 
nosity wave that triggers expansion of the outer layers, and then rapid contraction. 
Although some quantitative differences are apparent, significant differences in be- 
haviour do not occur until the end of the contraction phase. At this point, the further 
evolution of the protostar is governed by the size of the accretion rate. For accretion 
rates greater than some critical value M cr j t , the luminosity of the protostar becomes 
large enough to halt the accretion. On the other hand, for M < M cr jt, the lower accre- 
tion luminosity means that the total luminosity of the protostar remains below L^, 
an d accretion contin ues unabated. 

Omukai & Palla ( 2003} solve for M c rit by equating the total luminosity of a zero- 



age main sequence Pop. Ill protostar (including accretion luminosity) with the Ed- 
dington luminosity, and find that 

M crit ~4xlO- 3 M yr- 1 , (77) 

coincidentally close to Mm- m principle, one would expect M cr ; t to have a depen- 
dence on the current mass of the protostar, but in practice, |Omukai & Palla| {2003 ) 
show that this dependence is weak and may be neglected. 

Finally, Omukai & Palla ( 2003 ) show that in the time-dependent accretion model, 



the key factor is the size of the accretion rate at the end of the contraction phase. If 
this is greater than M cr ; t , then one would expect accretion to be halted, while if it is 



less than M C rit then accretion can continue. In practice, Omukai & Palla ( 2003 ) show 



that if one adopts the Abel, Bryan & Norman (2002) estimated accretion rate, then 



M < Merit, implying that accretion can continue even once the protostar reaches the 
zero-age main sequence. 

The main limitation of the approach outlined above is the neglect of the effects of 
rotation. In reality, rotation can have profound effects on stellar structure and evolu- 
tion, particularly for massive stars ( |Maeder & Meynet] |2000| l, and it will also have 
a large influence on how matter reaches the protostar in the first place. The first de- 
tailed study of the pre-main sequence evolution of a Pop. Ill protostar to account for 
the effects of rotation was carried out by Tan & McKee (2004). In contrast to pre- 



vious authors, they did not assume spherical symmetry. Instead, they assumed that 
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a protostellar accretion disk would form, and fixed the size of the disk by assum- 
ing angular momentum conservation within the supersonic portion of the accretion 
flow. They used the polytropic accretion flow model described in the previous sec- 
tion to compute the accretion flow onto the disk. To solve for the disk structure, 
they made use of the standard theory of steady, thin viscous accretion disks (as out- 
lined in Shakur a~& Sunyaev|1973| l, with a spatially constant viscosity parameter a. 



As sources for a, they considered the magnetorotational instability Balbus & Haw- 



ley (1991 1998] > and gravitational instability. With the disk structure in hand, they 



could then solve for the structure of the protostar itself, using a modified version 
of an approach developed by Nakano, Hasegawa & Norman (1995T) and Nakano 



|et al.| ( [2000l ). In the zero angular momentum case, |Tan & McKee| ( |2004"l ) show that 
they successfully reproduce the previous results of Stahler, Palla & Salpeter ( 1986a) 
and |Omukai & Palla| ( [200T1 12003) . In more realistic models, |Tan & McKee| ^2004 ) 
show that the presence of an accretion disk has little influence on the evolution of 
the protostar, which still evolves through the same progression of adiabatic growth, 
terminated by the emergence of a luminosity wave, followed by rapid contraction to 
the ZAMS. However, [Tan & McKee| ( |2004| l do find that the photosphere surround- 
ing the protostar behaves very differently in this case than in the spherical case. 
Because most of the gas accretes onto the protostar via the disk, the gas density is 
significantly reduced in the polar regions. Consequently, the optical depth of these 
regions is also significantly reduced, with the result that the flow becomes optically 
thin early in its evolution. For example, in the model with /kep = 0.5, the photo- 
sphere vanishes once the protostellar mass reaches 1 M and does not subsequently 
reappear. |Tan & Mc Kee (2004) argue that this may have a major influence on the 
effectiveness of radiative feedback from the protostar, a topic that we will return to 
in the next section. 



3.1.3 Feedback effects 

Accretion of gas onto the protostar liberates a significant amount of energy, with 
most of this energy being emitted from regions close to the protostellar surface. This 
can be shown very simply by considering how the gravitational potential energy of 
a test mass changes as we move it close to a protostar of mass M t and radius . 
At a distance of the gravitational potential energy of a fluid element with mass 
dM is 

W 

while at the protostellar surface it is 

W 

Therefore, the amount of energy that must be dissipated by the fluid element as it 
moves from 2R* to is as large as the amount that it must have dissipated while 



GM^dM 
2R* 



GM^dM 
R* 



(78) 



(79) 
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falling in from R ^> R* to 2R*, or in other words, half of the total binding energy 
dissipated by the gas is dissipated while its distance from the protostellar surface 
is less than R*. In addition, once the protostar reaches the main sequence, it will 
start generating additional energy in its own right, via nuclear fusion. The energy 
that is released in the vicinity of the protostar is therefore quite considerable, and 
it is reasonable to suppose that this will have some effect on the behaviour of the 
surrounding gas. It is therefore not surprising that considerable attention has been 
paid to the issue of protostellar feedback in the context of Pop. Ill star formation. 

In order for the protostar to substantially reduce the rate at which matter flows 
onto it, it must be able to transfer a significant amount of energy and/or momentum 
to the infalling gas. The various mechanisms by which this can be accomplished 
fall under two broad headings: mechanical feedback, where the protostar transfers 
energy and momentum to some form of outflow, which subsequently transfers it to 
the infalling material, and radiative feedback, where radiation from the protostar 
transfers energy and momentum directly to the infalling gas. 



Mechanical feedback 

In the local Universe, stellar winds are an almost ubiquitous phenomenon, and play 



an important role in the evolution of the most massive stars (Chiosi & Maeder 



1986). However, there are good reasons to expect that metal-free stars will be much 
less effective at driving winds than the roughly solar metallicity stars that we are 
familiar with in the Milky Way. Strong stellar winds are invariably radiation-driven, 
and at solar metallicities, the largest contribution to the radiative acceleration of 
the gas comes from the absorption and scattering of ultraviolet photons in the lines 



of the many metal atoms and ions present in the outflowing gas ( |Castor, A bbott, 



|& Klein| |1975| l. In metal-free gas, on the other hand, the only significant sources 
of opacity within an outflow will come from the lines of He + (atomic hydrogen 
is typically fully ionized), and from Thomson scattering by free electrons. These 
provide orders of magnitude less radiative acceleration per unit luminosity than do 
the metal lines in a solar metallicity gas, and hence one can show that a metal-free 
Population III star can produce a line-driven wind only if the stellar luminosity is 
already very close to the Eddington limit (Kudritzki 2002). 

Of course, as a Population III star evolves, it will not remain metal-free. It will 
start to produce carbon, nitrogen and oxygen internally once the stellar core begins 
to burn helium, and if the star is rotating, these elements can become well-mixed 
within the star (Meynet, Ekstrom & Maeder, 2006). This will provide an additional 
source of opacity in the stellar atmosphere which may allow the most massive Pop- 
ulation III stars to produce a weak CNO-driven wind ([Rrticka & Kubat[ 2009). 



However, the mass-loss rate will be small, and the fraction of the stellar mass that 
can lost in this way is unlikely to be larger than about 1%. 

It is also possible that very massive Population III stars with luminosities close 
to the Eddington luminosity may produce eruptive, continuum-driven winds, similar 
to those we see coming from nearby luminous blue variables (LBVs) such as rj Car 
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(Smit h & Owocki| [2006). However, as the cause of these LBV eruptions is not yet 
fully understood even for nearby objects, it is difficult to say with certainty whether 
they will actually be produced by Pop. Ill stars. More work on this topic is clearly 
necessary. 

Finally, mechanical feedback can also be generated in the form of hydrodynam- 
ical or magnetohydrodynamical jets or outflows. We have already discussed the 



magnetically-driven disk winds produced in the Machida et al. ( 2006 2008 



ulations, which are able to eject roughly 10% of the infalling gas from the disk. 
Although, as we noted previously, these simulations only modelled the very earliest 
stages in the formation of the protostellar accretion disk, their value for the mass 
ejection rate is in good agreement with the predictions of a semi-analytical study of 
Pop. Ill disks and outflows carried out by |Tan & Blackman| ( [20 04 ). If this value is 
correct, then it implies that the reduction in the protostellar accretion rate brought 
about by these outflows is small, and hence that they will not significantly limit the 
final stellar mass. However, one should bear in mind that their interaction with the 
star-forming halo on larger scales has not yet been modelled in any detail, and hence 
it is difficult to be certain regarding their final impact. 



Radiative feedback 

There are several different forms of radiative feedback that could potentially affect 
the accretion of gas by a Pop. Ill protostar. First, if the radiation is absorbed or 
scattered, then it will exert a force on the gas. If this force is comparable to or larger 
than the gravitational force acting on the gas, then it may suppress accretion onto 
the protostar, or even prevent it completely. Second, radiation may destroy the H2 
molecules responsible for cooling the gas. In the absence of cooling, the gas will 
evolve adiabatically, which again may reduce the rate at which it can be accreted. 
Third, the radiation may heat the gas. If radiative heating raises the gas temperature 
to a point at which the thermal energy of the gas exceeds the gravitational binding 
energy of the system, then this again will strongly suppress accretion. 

In local star-forming regions, the first of these three forms of radiative feedback 
is believed to be the most important. Radiation pressure exerted on infalling dust 
grains by radiation from the protostar results in a substantial momentum transfer 
to the dust, and from there to the gas, since the dust and gas are strongly coupled. 
In spherically symmetric models, the radiative force exerted by the radiation on 
the dust can be strong enough to bring accretion to a complete halt (Wol fire &] 
|CassinelEl |1987| ). In primordial gas, there is no dust, and so this process cannot 
operate. However, radiation pressure can also work directly on the gas, and so it 
is worthwhile investigating whether this process is likely to significantly suppress 
accretion. 

Let us start by assuming that the bolometric luminosity of the protostar is given 
by the Eddington luminosity 
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^Edd 



AkGM^c 



(80) 



where is the protostellar mass, and Kp = Opj '«p — 0.4 cm 2 g~ 1 is the opacity due 
to Thomson scattering for a fully ionized gas composed of pure hydrogen, with Op 
the Thomson scattering cross-section of the electron and ra p the mass of the proton. 
In this case, then we know from the definition of the Eddington luminosity that the 
radiative force exerted on a fluid element will be equal to the gravitational force 
exerted on it by the protostar when the opacity of the fluid element is equal to Kp. 
More generally, we can write the ratio of the forces acting on the fluid element as 



^rad 



^Edd KT ' 



(81) 



where F a d is the radiative force, f grav is the gravitational force, L t is the protostellar 
luminosity, and K is the mean opacity of the fluid element. Since the protostar is 
unlikely to be stable if > Lgdd, this implies that in order for the radiative force 
to significantly affect the gas, it must have a mean opacity K ~ Kp or higher. In 
practice, the luminosity of a Pop. Ill protostar before it reaches the main sequence 



will often be significantly less than the Eddington luminosity (see e.g. |Smith et al. 
201 l| l, in which case an even higher mean opacity is required. 



The mean opacity of metal-free gas has been computed by a number of authors, 
most recently by Mayer & Duschl (2005b). They present tabulated values for both 
the Rosseland mean opacity 



Jo(dB v /dT)K- 1 dv 
J~(dB v /dT)dv ' 



and the Planck mean opacity 



Kp 



IoB v dv ' 



(82) 



(83) 



where K v is the frequency-dependent opacity and B v is the Planck function. For our 
purposes, we are most interested in the Planck mean. Strictly speaking, this Planck 
mean opacity is the same as the mean opacity in Equation|ST|only if the protostar has 
a black-body radiation field and a photospheric temperature that is the same as the 
gas temperature, and in general this will not be the case. However, if the protostar 
is still in the pre-main sequence phase of its evolution, it will have a photospheric 
temperature T p ~ 6000 K ( Stahler, Palla & Salpeter 1986a) and a spectrum that does 



not differ too greatly from a black-body, while the temperature of the surrounding 



2011b 



gas will typically be of the order of 1000-2000 K or higher (Clark et al 
Greif et al. 2011a Smith et aLj 201 1) . In these conditions, the error we make by 



using the Planck mean opacity in Equation 81 should not be excessively large. 

Ther e are two regimes in which the tabulated values of Kp in Mayer & Duschl 
(2005b I exceed Kp. The first occurs at very high densities (n > 10 22 cm -3 ), where 
Kp > Kp for a wide range of temperatures. However, these extreme densities are 
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only reached within the protostar and hence this regime is of no relevance when 
we are considering feedback from the protostar on the surrounding gas. The second 
regime in which Kp grows to the required size is at temperatures above 8000 K, for 
a wide range of densities. At these temperatures, the dominant source of opacity is 
the scattering of photons in the Lyman series lines of hydrogen, primarily Lyman- 
a. The effects of Lyman-a radiation pressure in metal-free gas were considered by 
Oh & Haiman (2002 ), in the context of the formation of massive star-forming mini- 



halos with virial temperatures T > 10 K. They argued that the Lyman-a photons 
produced by the cooling of the hot gas would not be important (see also |Rees & 



Ostriker||1977j >, but that the Lyman-a photons produced by a massive star and its 



associated HII region would have a pronounced effect on the gas, and could signif- 
icantly delay or even halt the inflow of the gas. However, they did not carry out a 
full quantitative investigation of the effects of Lyman-a radiation pressure. More 
recently, this issue was revisited by McKee & Tan (2008 J, who studied it in some 
detail. They found that in a rotating flow, most of the Lyman-a photons would even- 
tually escape along the polar axis of the flow, as it is here that the optical depths are 
smallest. They showed that if the rotational speed of the gas were at least 10% of the 
Keplerian velocity, then Lyman-a radiation pressure would be able to reverse the di- 
rection of the flow along the polar axis once the protostellar mass reached 20 M . 
The radiation would therefore blow out a polar cavity, allowing more Lyman-a pho- 
tons to escape. This prevents the radiation pressure from rising further, and McKee 
& Tan argue that it never becomes large enough to significantly affect the inflow of 
gas from directions far away from the polar axis (e.g. from the accretion disk). For 
this reason, they conclude that Lyman-a radiation pressure is unlikely to be able to 
significantly reduce the protostellar accretion rate. 

Let us now turn our attention to the second form of radiative feedback men- 
tioned above: the photodissociation of H2 and the consequent dramatic reduction 
in the cooling rate. As we have already discussed, at early times the photospheric 
temperature of the protostar is too low for it to produce significant quantities of far- 
ultraviolet radiation, and hence radiation from the protostar does not significantly 
affect the H2. Once the protostar reaches the main sequence, however, it can become 
a significant source of far-ultraviolet radiation, provided that it has a mass greater 
than around 15 M Q ([McKee & Tan| |2"0"08"1 >. Studies by |Omukai & Nishi| ( [T999] > and 
|Glover & Brand| ( |200l| l considered the effect that this radiation would have on the 
H2 surrounding the protostar, and showed that the time required to photodissociate 
the H2 would be significantly less than the lifetime of the protostar. The removal of 
the H2 from the gas means that it is no longer able to cool effectively at temperatures 
T < 10 4 K, and hence one would expect that as the H2 in the accreting gas is de- 
stroyed, the gas will begin to evolve adiabatically until it reaches this temperature. 
McK ee & Tan| ( |2008| l consider whether this switch to adiabatic evolution is suffi- 
cient to halt accretion, and conclude that it is not. If no protostar were present, then 
the switch to adiabatic evolution would be enough to stabilize the gas and prevent 
further collapse. The presence of the protostar, however, serves to destabilize the 
gas, allowing accretion to continue even when the evolution of the gas is fully adi- 
abatic. iM^K^e&TSn] d2008 ) use the treatment of protostellar accretion introduced 
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in |Fatuzzo, Adams & Myers] ( |2004| l to investigate the issue numerically, and show 
that an increase in the effective adiabatic increase of the gas from y e g =1.1 (which 
approximately characterizes the temperature evolution of the gas at n > 10 4 cm~ 3 ; 



see e.g. Omukai & Nishi 1998| l to y e ff = 5/3 reduces the accretion rate by only 20%. 



The third possible form of radiative feedback involves the heating of the sur- 
rounding gas by radiation from the protostar. If the temperature of the gas can be 
increased to a point at which its thermal energy exceeds its gravitational binding en- 
ergy, then it will no longer be gravitationally bound to the protostar, and hence will 
not be accreted. A convenient way to quantify the relative importance of thermal 
and gravitational energy is to compare the sound-speed of the gas with the escape 
velocity of the system, v eS c: gas with c s > v eS c will not be gravitationally bound. 

For an isolated protostar of mass M„, we can write v esc at a distance R from the 
protostar as: 

(84) 

R 

where G is the gravitational constant. If we rewrite this expression in more conve- 
nient units, we find that 

v esc ~ 4.2 I ^] V2 ( ^ kms- 1 . (85) 



im s j v i00AU , 

For a primordial, fully molecular gas, c s = 4.2 km s -1 at a temperature T ~ 3400 K, 
and hence gas within 100 AU of a one solar mass protostar must be heated up to 
a temperature of thousands of Kelvin in order to unbind it. At larger distances, the 
required temperature would appear at first to be much smaller, but the reader should 
recall that this expression is for an isolated protostar, i.e. one which is not sur- 
rounded by gas. It is therefore only valid when the protostellar mass is much 
larger than the mass of gas within a distance R of the protostar, and once we start 
considering scales R ^> 100 AU, this is unlikely to be a good approximation. If we 



include the influence of this gas by replacing in Equation 84 by M tot = +M gas , 
and use the facts that prior to star formation, the mass enclosed within a sphere of 
radius 100 AU is roughly 5 M E , and increases at larger distances as M enc °c then 
at distances R > 100 AU, we have 



( R 



-o.i 



:9.4( — - — | kms" 1 . (86) 

ViooAuy v 

In other words, once we account for the mass of the infalling gas in addition to the 
mass of the protostar, we find that the escape velocity is of the order of lOkms -1 , 
with little dependence on the distance from the protostar. An escape velocity of this 
order of magnitude corresponds to a gas temperature of order 10 4 K. This imme- 
diately tells us that heating of the gas by radiation from the protostar during the 
pre-main sequence phase of its evolution is unlikely to significant affect the accre- 
tion rate due to the low photospheric temperature of the protostar during this phase 
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- clearly, a protostar with an effective temperature of 6000 K will not be able to heat 
up distant gas to a temperature of 10000 K. On the other hand, once the protostar 
reaches the main sequence, its photospheric temperature will sharply increase, and 
hence it may be able to heat up the surrounding gas to a much higher temperature. 
In particular, if the protostar is massive enough to emit a significant number of ion- 
izing photons while on the main sequence, then it will easily be able to produce 
temperatures in excess of 10 4 K within the gas that it ionizes. 

The idea that the formation of an HII region may strongly suppress or completely 
terminate protostellar accretion was discussed long ago in the context of present- 



day star formation (see e.g. Larson & Starrfield 1971 1, but has recently been re 



examined by several authors in the context of primordial star formation. On large 
scales (R > 0. 1 pc), the behaviour of an HII region produced by a Pop. Ill star is rel- 
atively simple. The radial density profile of the gas on these scales is approximately 
p «: R~ 2 2 , and hence the density falls off too quickly to trap the HII region within 
the minihalo (see e.g WhalenTAbel & Norman| [2004 |Alvarez, Bromm & Shapiro 



20061 |Abel, Wise & Bryanj |2007HYoshida et ai.||2007) . The ionization front there 



fore expands rapidly, as an R-type front, with a velocity that is controlled by the rate 
at which ionizing photons are being produced by the star. In addition, if we are con- 
sidering Pop. Ill star formation within one of the first star-forming minihalos, then it 
is easy to show that sound speed of the gas within the HII region will be higher than 
the escape velocity of the minihalo. Consequently, the ionized gas begins to flow out 
of these small minihalos, significantly reducing the mean gas density. It is therefore 
clear that once the HII region reaches a size of 0.1 pc or above, it will act to prevent 
any further infall of gas from these scales onto the protostar. However, this leaves 
unanswered the question of how long it takes for the HII region to expand to this 
scale. 



In the case of steady, spherically-symmetric infall, Omukai & Inutsuka (2002) 
showed that in order for the HII region to avoid being trapped on scales close to the 
protostar, the flux of ionizing photons must exceed a critical value 

/ \ — 1 / \ 2 

52 I R m \ I M * \ -1 



^ = 6 - 4xlo ~UofJ UomJ s < (87) 

where R m is the inner radius of the HII region, which we can take to be equal to the 
radius of the massive star. Given reasonable values for /?; n and M*, this expression 
yields a value for N CI it that is much larger than the number of photons that will 
actually be produced by any massive star, leading Omukai & Inutsuka (2002]) to 
conclude that the HII region would remain trapped close to the star. However, this 
conclusion depends crucially on the assumed spherical symmetry of the flow. In 
the more realistic case in which our protostar is surrounded by an accretion disk, 
McKee & Tan (2008]) show that the HII region can expand in all directions other 



than those close to the midplane of the disk once the stellar mass reaches a value 
of around 50— 100 M©, where the precise value required depends on how rapidly the 



gas is rotating. McKee & Tan (2008 1 also show that the accretion disk can survive 



for a considerable period after the HII region has broken out, and that the protostar 
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will stop accreting from the disk only once the rate at which gas is lost from the disk 
by photoevaporation exceeds the rate at which fresh gas is falling onto the disk. In 
their models, this occurs forM* ~ 140 M Q , leading them to conclude that radiative 
feedback from the protostar on the surrounding gas cannot prevent the protostellar 
mass from becoming very large. On the other hand, initial attempts to model this 
process in 2D or 3D ( |Hosokawa eTaL] |20~TTj |Stacy, Greif & Bromm| [2012] ) find 
significantly smaller final masses, M* ~ 30-40 M , although these detailed models 
have so far explored only a small part of the potential parameter space. 



3.2 The fragmentation model 
3.2.1 Early studies 

The first simulations of primordial gas to make use of sink particles were the SPH 
simulations of Bromm, Copp i & Larson| ([T999 2002 ). They studied the formation 
of isolated dark matter minihalos and the cooling and gravitational collapse of gas 
within them using a somewhat idealized set of initial conditions. At an initial red- 
shift z = 100, they created a spherical, uniform density region containing both gas 
and dark matter, and with an initial velocity field corresponding to the Hubble ex- 
pansion. The density of this spherical region was taken to be higher than the cosmo- 
logical background density, and was fixed such that the region would gravitationally 
collapse and virialize at a specified redshift, chosen to be z = 30 in most of the mod- 
els that they examined. Small-scale structure was introduced into the dark matter 
distribution by perturbing the particles slightly from their initial positions using the 



Zel'Dovich ( 1970 1 approximation. The amplitudes of these random perturbations 
were fixed such that the small-scale density structure in the dark matter would begin 
to evolve in the non-linear regime at the virialization redshift. Both the dark matter 
and the gas were also assumed to be in solid body rotation, with some specified 
angular velocity. 



Bromm, Coppi & Larson| examined several different choices for the halo mass 



and initial angular velocity of the gas, and showed that starting from these initial 
conditions, the gas and dark matter would initially collapse in a similar fashion, 
but that the gas would subsequently form H2, dissipate energy, and sink to the cen- 
ter of the minihalo. In most of the cases they studied, the gas would then form a 
rotationally-supported disk, with a radius of order 10 pc. This disk would then break 
up into clumps with masses M c \ ~ 100-1000 M , comparable to the Jeans mass in 
the disk. As these clumps were gravitationally unstable, they of course underwent 
gravitational collapse, and |Bromm, Coppi & Larson ( 2002[ ) therefore introduced 



sink particles to represent clumps that collapsed to densities greater than 10 8 cm -3 
in order to avoid the timestep constraints discussed earlier, allowing the further evo- 
lution of the clumps to be studied. 

Unfortunately, the fragmentation observed by Bromm, Cop pi & Larson| in their 
simulations is probably not realistic. One major problem lies in their choice of initial 



50 



Simon Glover 



conditions, specifically in their use of solid-body rotation. Although the total angu- 
lar momentum of the gas and dark matter in their simulations is comparable that 
measured for minihalos in more realistic cosmological simulations (|Jang-Condell 



&Hernquist 2001 Davis & Nataraj an 2010[), their adoption of solid-body rotation 



leads to the gas having an incorrect radial profile for this angular momentum. This 
causes the collapse of the gas to be considerably more ordered than it would be in 
a real minihalo, leading to the formation of an over-large disk. Disks of this kind 
do not appear to form in simulations of small star-forming minihalos that start from 



more realistic cosmological initial conditions (e.g. Abel, Bryan & Norman 2002 



Yoshida et al. 2006 ). A second major problem lies in the neglect of stellar feedback. 



Within the disk, the dynamical timescale is of the order of a million years, which is 
much longer than is needed for a massive Pop. Ill star to reach the main sequence. 
Therefore, if a massive star forms within the first clump to be produced within the 
disk, the radiation from this star may well be able to photodissociate the H2 in the 



1999, Glover & Brand 



disk before a second clump can form (Omukai & Nishi 
12001) . 

The next attempt to use sink particles to study the formation of Pop. Ill stars was 
made by Bromm & Loeb ( 2004)1. The init ial setup of their simulation was similar 
to that used by Bromm, Coppi & Larson (2002), but to gain improved resolution 



in the centre of the minihalo, they used a technique called particle splitting (Kit- 



|sionas &~W hitworth 2002; Bromm & Loeb||2003| l. The evolution of the minihalo 
was followed until cold, dense gas started to build up in the centre of the halo. 
The simulation was then paused, and the gas within a radius of 3.1 pc of the cen- 
tre of the minihalo (corresponding to roughly 3000 M of material) was resam- 
pled using SPH particles with much smaller masses, using the resampling tech- 
nique described in Bromm & Loeb (2003 1. The mass resolution within this resam- 
pled region was thereby improved from M res = 200 M to M res = 4 M . Bromm & 
Loeb then restarted the simulation, and followed the further gravitational collapse 
of the gas within this central, higher resolution region until the gas density reached 
n ~ 10 12 cm~ 3 , at which point a sink particle was created. They then followed the 
accretion of gas onto this sink for roughly 10 4 years, as we have already described 
above. Bromir f& Loeb| ( 2004| l found no evidence for fragmentation within the cen- 
tral clump of dense gas, but did note that a second dense clump formed nearby, with 
a final separation from the star-forming clump of roughly 0.25 pc. However, the 
free-fall collapse time of this clump was about 3 Myr, and so it was unclear whether 
it would survive for long enough to form a second star, or whether it would be 
destroyed by negative feedback from a massive star forming within the first clump. 



3.2.2 The importance of turbulence 

Although the |Bromm & Loeb (2004i study undoubtedly represented a significant 
step forwards in resolution compared to |Bro mm, Coppi & Larson (2002 ), it still had 
a mass resolution which was more than two orders of magnitude greater than the ac- 
tual size of a Pop. Ill protostar at the moment that it forms, and hence it was unable 
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to investigate the behaviour of the gas on scales smaller than about 100 AU. The first 



work using sink particles that did manage to probe this regime was Clark, Glover 



|& Kle ssen (2008 ). Although the main focus of their study was on the fragmentation 



brought about by dust cooling in low metallicity systems (see e.g. Omukai et al. 



|2005 1 [Schneider et al.|2006 or Dopcke et al.|2011 for more on this topic), they also 
studied the behaviour of the gas in the Z = case for the purpose of comparing it 
with the results on their low metallicity runs. As the initial conditions for their sim- 
ulations, Clark, Glover & Klessen ( 2008| > considered a uniform density cloud, with 
a mass of 500 M Q , a radius of 0.17 pc and a number density of 5 x 10 5 cm~ 3 . The 
gas within this cloud was given a low level of initial turbulence, with a turbulent en- 
ergy equal to 10% of the gravitational potential energy, and was also assumed to be 
rotating uniformly, with an initial rotational energy equal to 2% of the gravitational 
potential energy. Two different simulations were performed, with different numbers 
of particles: a low resolution simulation that used only two million SPH particles, 
and hence had a mass resolution of 0.025 M , and a high resolution simulation that 
used 25 million particles, corresponding to a mass resolution of 2 x 10~ 3 M . Aside 
from the somewhat artificial initial conditions, the main simplification made in these 
simulations was the use of a tabulated equation of state to follow the thermal evolu- 
tion of the gas. The results of the Omukai et al. ( 2005 ) one-zone model were used to 
derive internal energy densities and thermal pressures for the gas at a range of dif- 
ferent densities, and this data was then used to construct a look-up table that could 
be used by the SPH code to compute the internal energy and pressure corresponding 
to a given gas density. 

Clark et al. followed the collapse of the gas in their simulation down to a physical 
scale of less than an AU (corresponding to a gas density of over 10 16 cm -3 ). Regions 
collapsing to even smaller scales were replaced by sink particles, created using the 
standard [Bate, Bonnell & Price ( 1995) prescription. Clark et al. showed that at the 
point at which the first sink particle formed, the radial profiles of quantities such as 
the gas density or the specific angular momentum were very similar to those found 
in previous studies of Pop. Ill star formation that were initialized on cosmological 
scales (e.g. |Abel, Bryan & Norman] [2002[ |Yoshida et ah] [2006] ). They noted that 
at this point in the simulation, there is no sign of any fragmentation occurring, and 
argued that if the simulation were stopped at this point (as would be necessary if the 
sink particle technique were not being used), one would probably conclude that the 
gas would not fragment, but would merely be accreted by the protostar. However, 
they show that this is not what actually happens when the simulation is continued. 
Instead, the gas fragments, forming 25 separate protostars after only a few hundred 
years. Clark et al. stopped their simulations after 19 M of gas had been incorpo- 
rated into sink particles, and showed that at this point the protostars have masses 
ranging from 0.02 M to 5M , but that the mass distribution is relatively flat, with 
most of the mass locked up in the few most massive protostars. There is no signif- 
icant difference between the mass function of sinks in the low and high resolution 
calculations, suggesting that fragmentation is well-resolved in both cases. 

This is an intriguing result, but several reasonable concerns could be raised re- 
garding the numerical technique adopted by Clark, Glover & Klessen ( 2008 ). First, 
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the initial conditions for the gas are an idealized version of what one would find 
within a real star-forming minihalo, and although there are indications that the 
gas loses its memory of the initial conditions prior to fragmentation occurring, in- 
evitably a few doubts remain. Second, and more importantly, the use of a tabulated 
equation of state represents a major simplification of the thermal evolution of the 
gas, and one which may make fragmentation more likely to occur. For example, this 
technique does not allow one to model the formation of the hot, shocked regions 
noted by Turk, Norman & Abel (20101 in which much or all of the H2 is dissoci- 
ated, and it is likely to underestimate the temperature of gas falling in at later times, 
when the typical infall velocity is larger than during the initial assembly of the pro- 
tostar. The Clark et al. calculation also neglects the effects of radiative feedback 
from the protostars, and assumes that protostars do not merge, even if they come 
within sub-AU distances of each other. 

In a follow-up study, Clark et al. ( |201 la| > addressed one of these concerns - 
the use of a tabulated equation of state - by performing simulations that replaced 
this with a detailed treatment of the chemistry and cooling of primordial gas. In 
their study, they investigated the role that low Mach number turbulence might play 
in triggering fragmentation in the gas by performing a set of simulations of the 
collapse of unstable Bonnor-Ebert spheres. They considered three initial configura- 
tions: a 1000 Mq cloud with an initial temperature of 300 K; a 150 Mg cloud with 
an initial temperature of 75 K; and a 1OOOM cloud with an initial temperature of 
75 K. In each case, the central density of the Bonnor-Ebert sphere was taken to be 
n c = 10 s cm~ 3 . The first set of initial conditions were intended to correspond to the 
conditions that one would expect to find within one of the first star-forming mini- 
halos, while the second set were intended to correspond to the conditions within a 
minihalo dominated by HD cooling. Simulations with the third set of initial con- 
ditions were run to allow the effects of lowering the temperature and lowering the 
total mass to be distinguished. Within the Bonnor-Ebert spheres, a turbulent velocity 
field was imposed, with a three-dimensional RMS velocity Ziv tur t,. 

Clark et al. ( 201 la\ did not claim that this was a completely accurate model of 



the physical state of the gas within a real star-forming minihalo. Instead, they treated 
this study as a kind of physics experiment, allowing them to investigate the effect 
of varying a single important parameter - the turbulent kinetic energy - without 
varying any of the other parameters in the problem, something that it would not 
be possible to do if using initial conditions derived from a cosmological simula- 
tion. For the first two setups described above, they performed four simulations, with 
^ y tm-b = 0.1,0.2,0.4 and 0.8c s , respectively, where c s was the initial sound speed. 
For the third setup (the large, low temperature clouds), they performed only two 
simulations, with 4v tur b = 0.4 and 0.8 c s , respectively. The clouds were modelled 
using two million SPH particles in each case, yielding a mass resolution of 0.05 M 
for the 1000 M e clouds and 0.0075 M for the 150 M clouds. Sink particles were 
created once the gas density exceeded 10 13 cm~ 3 , and the sink accretion radius was 
20 AU. 

|Clark et aL"1 ( |201 la| found that fragmentation occurred in almost all of their sim- 
ulated clouds. In the case of the massive, warm clouds, the only case in which frag- 
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mentation did not occur was the simulation with 4v tur b = 0.1 c s . In this simulation, 
the gas simply collapsed to form a single, massive protostar. In the simulations with 
larger turbulent energies, however, the formation of the first protostar was followed 
within a couple of hundred years by the fragmentation of the infalling gas and the 
formation of a significantly larger number of protostars. The relationship between 
the turbulent energy and the degree of fragmentation is not straightforward: the 
^Vturb = 0.4 c s run fragmented more than the 4v tur b = 0.2 c s , as one might expect, 
but the 4v tur b = 0.8c s run fragmented less than the Avturb = 0.4c s run (although still 
more than the 4v tur b = 0.2 c s run). Clark et al. hypothesize that this difference in be- 
haviour is due to the amount of angular momentum retained within the collapsing 
region, which in this case was larger in the Avturb = 0.4 c s run than in the other runs, 
but note that this may not always be the case, and that a much larger series of realiza- 
tions of the turbulent velocity field would be needed to make a definitive statement 
about the relationship between the turbulent energy and the degree of fragmentation 
(c.f. Goodwin, Whitworth, & Ward-Thompson 2004 who come to a similar con- 
clusion regarding present-day star formation). The total amount of mass accreted by 
the sinks is very similar in all four runs, and hence the runs that fragment more tend 
to form lower mass objects than the runs that fragment less. 

In the low-mass, colder clouds, Clark et al. find a much lower degree of frag- 
mentation, despite the fact that the initial number of Jeans masses in these clouds is 
the same as in the 1000 M , T — 300 K clouds. In this case, fragmentation occurs 
only in the Avturb = 0.2 c s and 4v tm -b = 0.8c s simulations, and only a small number 
of fragments are formed in each case. Clark et al. investigate whether this is due to 
the lower cloud mass by modelling the collapse of 1OOOM clouds with the same 
lower initial temperature, and find that although more fragmentation occurs in this 
case, the gas still fragments less than in the 1000 M , T = 300 K case. They sug- 
gest that this somewhat counterintuitive behaviour is due to the greater stiffness of 
the effective equation of state in the colder clouds. In both cases, the gas must heat 
up from its initial temperature at 10 5 cm~ 3 to a temperature of roughly 1000 K at 
10 10 cm~ 3 , and so when the initial temperature is lower, the gas must heat up more 
rapidly with increasing density, meaning that it has a larger effective adiabatic index. 
This makes it more difficult to generate the small-scale non-linear structures that are 
the seeds for later fragmentation, and also delays the collapse, allowing more of the 
turbulent energy to dissipate. A similar effect has previously been noted by Yoshida, 



Omukai & Hernquist (2007]l and Tsuribe & Omukai (2008), and calls into question 



the common wisdom that minihalos in which the cooling becomes HD-dominated 
will inevitably form lower mass stars. 



3.2.3 Models using cosmological initial conditions 



In order to establish whether the fragmentation seen in the Clark, Glover & Klessen 



( 2008 ) model was simply a consequence of the highly idealized initial conditions 



used in that study, several recent follow-up studies have re-examined the issue us- 
ing simulations initialized on cosmological scales (i.e. scales significant larger than 



54 



Simon Glover 



the virial radius of the minihalo). One of the first of these studies was carried out 



by Stacy, Greif & Bromm (2010). They first performed a medium resolution cos- 
mological simulation, which allowed them to determine the formation site of the 
first minihalo large enough to cool effectively and form stars. They then used a hi- 
erarchical zoom-in procedure (Navarro & White 1994| Tormen, Bouchet & White 



|1997| |Gao et ah] |2005[ ) to improve the resolution within a region centered on this 
formation site, allowing them to achieve a mass resolution of 1 .5 M R within th e 



Clark, Glover & Klessen 



(2008) 



centre of the star-forming minihalo^jln contrast to 
the thermal and chemical evolution of the gas was followed in detail during the col- 
lapse. Once the gravitationally collapsing gas reached a density of 10 12 cm~ 3 , it was 
converted into a sink particle, along with all of the gas within an accretion radius 
face = 50 AU. Stacy et al. show that following the formation of this first sink, the 
infalling gas collapses into a flattened disk. At a time t = 250 yr after the formation 
of the first sink particle, this disk has a radius of 200 AU, but it grows with time 
and has reached a radius of 2000 AU by t = 5000 yr, the end of the simulation. H2 
cooling allows the gas within the disk to remain at a temperature of roughly 1000 K, 
and the disk soon becomes gravitationally unstable, forming a second sink particle 
after roughly 300 yr, and a further three sinks after 4000-5000 years of evolution. 
At the end of the simulation, the first two sinks to form have become very mas- 
sive, with masses of 43 Mq and 13 respectively, while the three newer sinks still 
have masses ~ 1 M El , close to the resolution limit of the simulation. Stacy, Greif & 



|Bromm| ( |2010| ) do not include the effects of accretion luminosity in their simulation 
directly, but do assess its effects during a post-processing stage. They investigate 
the possible effects of radiation pressure, but show that this remains unimportant 
within their simulation throughout the period that they simulate, in agreement with 
our analysis above. 

The main drawback of the Stacy, Greif & Bromm (2010) study is their choice of 
mass resolution. At the point at which it fragments, the protostellar accretion disk in 
their simulation has a mass of roughly 35 M , and hence is resolved with only a few 
thousand SPH particles. This is two orders of magnitude smaller than the number 
of particles typically used to model gravitationally unstable accretion disks in the 



context of present-day star formation (see e.g. |Rice, Lodato & Armitage 2005 ) 



and it is questionable whether a few thousand particles is enough to properly model 



the dynamics of the disk. It is therefore possible that the results of Stacy, Greif & 



Bromm ( 2010| l may have suffered from some degree of artificial fragmentation. 



More recently, a study by Clark et al. ( 201 lb| > has dramatically improved the 
mass resolution used to model the build-up of a protostellar accretion disk around 
the first Population III protostar. Clark et al. use a similar basic strategy to Stacy 
et al., starting with a medium resolution cosmological simulation to identify the 
first star-forming minihalo, and then using a hierarchical zoom-in procedure to im- 



The value quoted here for the mass resolution of the|stacy, Greif & Bromm 



(2010 



simulation 



assumes that 100 or more SPH particles are required to resolve gravitationally bound structures, 
which is the typical resolution limit adopted in studies of present-day star formation. Stacy, Greif &| 
IBrommlpnO) assume that only 48 SPH particles are required, and hence quote a mass resolution 
that is roughly a factor of two smaller. 
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prove the resolution within the gas forming this minihalo. They run this zoomed-in 
simulation until the maximum density of the gravitationally collapsing gas reaches 
10 6 cm . At this point, they extract a spherical region containing 1OOOM of gas 
from the centre of the minihalo, and resimulate this region at much higher resolu- 



tion, using several nest ed levels of particle splitting fKitsionas & Whitworth 2002 



Bromm & Loeb 2003 1. At the final level of splitting, the particle mass is 1O- 5 M 
and the mass resolution is 10~ 3 M Q , several orders of magnitude better than in the 



Stacy, Greif & Bromm (2010 1 simulation. 

In addition to the extremely high mass resolution, the other main improvement 
in the Clark et al. study compared to previous work is its inclusion of the effects 
of accretion luminosity feedback directly within the simulation. To model this, the 
authors start by writing the bolometric accretion luminosity produced by a given 
protostar as 

^acc = — - — , (88) 

where M is the accretion rate onto that protostar, is the protostellar mass, and 
R* is the protostellar radius. As Clark et al. attempt to model only the first few hun- 
dred years of the evolution of the gas after the formation of the first protostar, i.e. 
a timescale much less than the protostellar Kelvin-Helmholtz relaxation timescale, 
they assume that the protostars remain in the adiabatic accretion phase of their evo- 



lution, with masses and radii that are related by (Stahler, Palla & Salpeter 1986a I 



M,\ a2 V M ^ 041 



R* = 26R & — = r ■ (89) 

°\M Q J 1.10 Wl.yr ' ) 

The only remaining uncertainty is then M, which can be directly measured within 
the simulation. Clark et al. next assume that the gas is heated by the accretion lumi- 
nosity at a rate 

where p is the mass density, r is the distance to the protostar, and Kp is the Planck 



mean opacity of the gas, calculated using the tabulated values given in Mayer & 
Duschl ( 2005b ). This expression assumes that the gas is optically thin, and hence 
will tend to overestimate the heating rate. 

Clark et al. model protostar formation using sink particles, which are created 
using the standard |Bate, Bonnell & Price ( 1995) algorithm, with a density threshold 



n t h = 10 cm . The sink accretion radius was set to 1.5 AU. At the point at which 
the first sink particle forms, the state of the gas in the Clark et al. simulation (e.g. the 
density profile and the distribution of specific angular momentum) is very similar to 
that seen in other high resolution simulations of Pop. Ill star formation. However, 
the authors show that at later times, a protostellar accretion disk begins to build up 
around the central protostar, just as in the Stacy et al. study. The significantly higher 
resolution of the Clark et al. simulation allows them to model the build up of this 
disk on scales much closer to the central protostar, and to resolve the disk with a 
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far larger number of SPH particles. The growth of the accretion disk is followed for 
around 100 years after the formation of the first protostar, and Clark et al. show that 
after around 90 years (corresponding to around 1 .5 orbital periods for the disk), the 
accretion disk begins to fragment, forming several low-mass protostars. At the time 
at which it fragments, the disk contains several solar masses of gas (and hence is 
resolved with several hundred thousand SPH particles), compared to around O.4M 
in the central protostar, and the disk radius is a few tens of AU. The state of the disk 
at the onset of fragmentation is illustrated in Figure [4] 




Fig. 4 The volume density of the protostellar accretion disk in the Clark et al. 201 la simula- 
tion immediately prior to fragmentation. The 'hole' in the centre of the disk corresponds to the 
location of the sink particle representing the first protostar to form, and occurs because we have 
not accounted for the mass in the sink particle when calculating the density. The accretion disk is 
gravitationally unstable, and has formed several spiral arms, one of which has begun to fragment. 



Clark et al. argue that the reason that the protostellar accretion disk fragments is 
that it is unable to transfer gas onto the protostar fast enough to keep up with the 
rate at which fresh gas is falling onto the disk. This causes the surface density of 
the disk to increase, which eventually results in it becoming gravitationally unstable 
and fragmenting. This argument can be made somewhat more quantitative if one 
treats the disk using the standard Sha kura & Sunyaev| ( |1973| l thin disk model, and 
hence writes the mass flow rate through the disk at a radius r as 
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M{r) = 3zac s (r)E(r)H(r), 



(91) 



where a is the viscosity parameter, and c s (r), H(r) and E(r) are the sound speed, 
disk thickness and surface density, respectively, at a radius r. Clark et al. use the 
values of c s , H and E provided by their simulation to show that M(r) is smaller than 
the accretion rate onto the disk for a wide range of radii, even if one adopts a = 1 
(which already presupposes that the disk is gravitationally unstable). The growth 
of the disk therefore appears to be unavoidable, leading to its fragmentation once 
portions of it develop a Toomre stability parameter Q < 1, where Q = c s K/nGE, 
with K here being the epicyclic frequency. 

Previous semi-analytical studies of the structure of Pop. Ill accretion disks came 
to a somewhat different conclusion regarding their stability, predicting that Q^S> 1 



throughout the disk ( Tan & McKee 2004 Tan & Blackman 2004 Mayer & Duschl 



2005a). However, these models neglected the effect of H2 cooling, motivated by the 
assumption that the H2 content of a Pop. Ill protostellar accretion disk is negligible 
(J. Tan, private communication). This assumption leaves H as the primary source 
of opacity at temperatures T ~ 7000 K and below (|Lenzuni, Chernoff & Salpeter 



1991 Mayer & Duschl 2005b). For gas in chemical equilibrium, the opacity of 



H decreases sharply with decreasing temperature, and consequently these models 
find the equilibrium temperature of the accretion disk to be high, T ~ 6000 K. In 
comparison, Clark et al. show that when the effects of H2 are taken into account, the 
characteristic temperature of the gas in the disk lies in the range of 1500 to 2000 K. 
The disks in these previous semi-analytical models could therefore transfer mass 
onto the protostar more rapidly than the disk in the Clark et al. simulation, owing to 
their higher sound-speed and larger thickness, but at the same time were also more 
stable against gravitational fragmentation. It is therefore not surprising that these 
previous studies predicted that the accretion disk should be stable, and highlights 
the crucial role played by H2 cooling in enabling disk fragmentation. 

In addition to the simulation described above, in which the value of M used 
to calculate the accretion luminosity was measured directly, Clark et al. also per- 
formed two additional simulations in which M was kept fixed, allowing them 
to investigate the role played by accretion luminosity heating. They considered 
cases with M = 10~ 3 M Q yr _1 (somewhat smaller than the measured value) and 
M = 10~ 2 M© yr _1 (larger than the measured value), and showed that as M (and 
hence L a cc) increase, the disk becomes warmer and thicker and takes longer to 
fragment. However, fragmentation still occurs in every case, demonstrating that 
accretion luminosity heating is unable to prevent the disk from fragmenting (c.f. 



Krumholz 2006 who argues that it plays a crucial role in suppressing fragmenta- 



tion in local star-forming systems). 



One drawback of the Clark et al. (2011b! study is that they examined only a 



single star-forming minihalo, and although they showed that the properties of this 
halo (e.g. mass, spin parameter, formation redshift) were similar to those of the 
minihalos modelled in previous studies of Pop. Ill star formation, nevertheless the 
suspicion remains that perhaps this particular minihalo was unusual in some way. 
This concern was addressed by |Greif et al. ( 201 la| >. They used the new moving- 
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mesh code AREPO ( Springe!, |2010| l to study Pop. Ill star formation in five different 
minihalos, using a sink particle algorithm to allow them to follow the evolution of 
the gas past the point at which the first protostar formed. In all five of the systems 
that they modelled, they found similar behaviour to that in the Clark et al. study: an 
accretion disk built up around the first protostar, became gravitationally unstable, 
and began to fragment after only a few years. These results suggest that Clark et al. 
were right to claim that disk fragmentation, and the resulting formation of Pop. Ill 
binary systems, or higher order multiple systems, is an almost inevitable outcome 
of Population III star formation. 



Greif et al. (201 la\ also examined the issue of whether the objects represented by 



the individual sink particles would truly survive as separate protostars, or whether 
they would simply merge into a single massive protostar as the system evolved fur- 
ther. They considered two different schemes for merging sink particles. In the stan- 
dard scheme, sink particles coming within a distance of 100 R© of each other were 
merged to form a single sink, provided that the total energy of the two-body system 
was negative. In an alternative model, utilizing what Greif et al. dub as "adhesive" 
sinks, the energy check was omitted, and sinks were always merged when within 
a distance of 100 R of each other. Greif et al. justify their choice of this critical 
distance in two ways: first, it is also the accretion radius adopted for their sinks, 
meaning that the gas flow on smaller scales close to the sinks is not resolved; and 
second, it is roughly equal to the maximum size of a pre-main sequence protostar 
predicted by the models discussed in Section [3.1.2| 

The majority of the protostars formed in the Greif et al. simulations have at least 
one close encounter with another protostar, but when the standard merger algorithm 
is used, many of these encounters result in a purely dynamical interaction, as the 
total energy of the protostellar pair is too large to allow them to merge. On the other 
hand, when the adhesive sinks are used, many of these encounters lead to mergers. 
Greif et al. show that although the total mass incorporated into protostars is roughly 
the same in both cases, the number of protostars that survive as individual objects is 
reduced by a factor of up to a few, and the mean protostellar mass is consequently 
higher. In both cases, the protostars have a relatively flat mass distribution, with 
most of the protostellar mass being accounted for by a small number of high-mass 
protostars. The protostars have a broad distribution of radial velocities, ranging from 
v ra( j ~ 1 kms -1 to v ra d ~ lOOkrns -1 , and in many cases the radial velocity is greater 
than the escape velocity of the central region of the minihalo. It is likely that this 
leads to a significant fraction of the Pop. Ill protostars escaping from the minihalo 
entirely, although Greif et al. do not follow their evolution for long enough to con- 
firm this. It is possible that these protostars will accrete very little additional gas 
once they escape from the high density region at the centre of the minihalo (see e.g. 
Johnson & Khochfar |201 1) , in which case their final masses would be very simi- 



lar to the masses that they have at the point at which they are ejected. Greif et al. 
show that in the standard case, a considerable number of these ejected protostars 
have masses M < 1 M©. If these protostars do indeed avoid accreting further gas af- 
ter their ejection, then they would have lifetimes that are comparable to the current 
age of the Universe. This suggests that it may be possible for some Population III 
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stars to survive until the present day. However, when the adhesive sinks are used, the 
number of ejected protostars with subsolar masses is greatly reduced, demonstrating 
that this conclusion is highly sensitive to our treatment of protostellar mergers. 

The majority of the Greif et al. simulations did not include the effects of the 
accretion luminosity generated by the collection of protostars. However, they did 
consider one case in which this was included, using the Clark et al. treatment with a 
fixed value for the accretion rate used in the determination of the accretion luminos- 
ity, M = 0.1 M yr~ 1 . The effect of this was to puff up the disk, causing fragmenta- 
tion to occur at a slightly larger distance from the initial protostar. However, despite 
the unrealistically high value adopted for M, fragmentation still occurred and the 
number of protostars that formed was barely affected. 

A more detailed study of the effects of accretion luminosity heating was carried 
out by |Smith et"aL] ( |2011[l. They use d Gadget to resimulate the central 2 pc of the 



minihalos simulated by Greif et al. (201 la i, starting at a time prior to protostar 
formation at which the peak density of the gas was around 10 9 cm~ 3 . Smith et al. 



(201 1 ) evolved these systems past the time at which the first protostar formed, and 



used sink particles with large accretion radii (r acc = 20 AU) to allow them to follow 
the dynamical evolution of the system for an extended period. For each of the five 
minihalos, they performed simulations both with and without accretion luminosity 
heating. When the accretion luminosity heating was included, it was treated in the 



same fashion as in Clark et al. ( 201 lb I. They found that in general, the effect of the 
accretion luminosity heating was to delay fragmentation and reduce the number of 
fragments formed. However, they also showed that the effect was relatively small, 
and had less influence on the number of fragments formed than did the intrinsic 
variation in halo properties arising from their different assembly histories. 



3.2.4 Open questions 

As the discussion in the previous section has shown, the past couple of years has 
seen a large increase in the number of simulations of Pop. Ill star formation that 
show evidence for fragmentation, suggesting that the older picture that had Pop. Ill 
stars forming in isolation with masses of 1OOM or more is in need of some revision. 
However, many aspects of the fragmentation scenario remain unclear. Some of the 
most important open questions are summarized below. 

• Are our treatments of optically thick H2 cooling and accretion luminosity heating 
adequate? 

The fragmentation of the gas observed in these simulations typically occurs at densi- 
ties at which H2 line cooling is optically thick, and hence may depend on the method 
used to account for the reduction that this causes in the cooling rate. At present, both 
of the methods in common usage represent relatively crude approximations, and it 
remains to be seen whether the behaviour of the gas will remain the same if a more 
accurate treatment is used. Similarly, the method currently used to treat the effects 
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of accretion luminosity heating also makes a number of major simplifications that 
may influence the outcome of the simulations. 

• What role do magnetic fields play? 

If a non-negligible magnetic field can be generated by the turbulent dynamo dur- 
ing the gravitational collapse of the gas, as discussed in Section [23] then this may 
influence the evolution and stability of the disk. The presence of a magnetic field 
may make the disk unstable via the magnetorotational instability ( |Tan & B lackman 



2004 Silk & Langer 2006 1, although the resulting mass transfer onto the proto- 
star is unlikely to be fast enough to prevent the disk from becoming gravitationally 
unstable. A more important effect may be magnetic braking of the infalling gas, 
which could act to significantly reduce the angular momentum of the gas reaching 
the disk (see e.g. Hennebelle & Ciardi 2009). Whether either of these effects can 



significantly suppress fragmentation remains to be determined. 
• How often do Population III protostars merge? 

Current simulations either ignore mergers entirely (e.g. |Clark et aL] |201 lb| [Smith 



et al. 201 l| l, or treat them using very simple approximations that do not properly 



account for the effects of tidal forces (e.g. |Greif et al.||2011a| l. However, it is clear 
from the results of the Greif et al. study that the method used to treat mergers has a 
significant influence on the number of protostars that survive, their mass distribution 
and their kinematics. Improving the accuracy with which protostellar mergers are 
treated within this kind of simulation is therefore an important priority. 

• Can we find some way to do without sink particles? 

The concerns outlined above regarding the way in which protostellar mergers are 
treated would be greatly ameliorated if we were able to run the simulations with- 
out sink particles, as in this case we would be able to model directly how the gas 
behaves on scales of the order of 100R©. To do this, however, it will be necessary 
to devise some scheme for treating these very small scales that does not fall foul of 
the Courant time constraint discussed previously. Greif et al. ( 2012| l have recently 



published the results of an initial effort along these lines, but were only able to 
simulate the first ten years of the evolution of the disk, owing the extremely high 
computational cost of the calculation. 

• How rapidly does H? photodissociation occur? 

Fragmentation is dependent on the cooling provided by H2 and does not occur in 



models of protostellar accretion disks that omit this effect (Tan & McKee 2004 



Mayer & Duschl 2005a 1. It is therefore highly probable that fragmentation will 
cease once the H2 has been photodissociated by Lyman- Werner band photons emit- 
ted from any massive protostars that form. What is not yet clear is how quickly 
this will occur. Mc Kee & Tan| ( |2008] l show that the number of Lyman-Werner band 
photons produced by a zero-age main sequence Population III star increases sharply 
with increasing stellar mass, before levelling off at a value S] w ~ 10 49 photons s _1 
forM* ~ 30 M e . If we assume that all of these photons are absorbed by H2 and that 
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20% of these absorptions lead to dissociation (Drain e~& Bertoldi| [l996), then the 
radiation from the star will photodissociate H2 at a rate = 0.1 M yr _1 , leading 
to complete removal of the H2 within only a few hundred years. However, it is likely 
that many of the available photons will not be absorbed by H2, either because they 
never coincide with one of the Lyman-Werner band lines, or because they escape 
along a direction in which most of the H2 has already been dissociated, or because 
they are absorbed by atomic hydrogen (Glover & Brand 2001 1, and so the disso- 



ciation time for the H2 could be significantly longer than suggested by this simple 
estimate, particularly once one accounts for the effects of three-body H2 formation. 

• How rapidly is the gas ionized? Does this completely suppress accretion, or sim- 
ply suppress fragmentation? 



As we have already discussed in Section 3.1.3 the most plausible mechanism for 



shutting off the supply of cold gas at the centre of the minihalo is the formation of an 
HII region whose thermal pressure is sufficient to expel most of the gas. However, 
only a few studies have looked at the interaction between the growth of the HII 
region and the evolution of the protostellar accretion disk. In particular, the issue 
has not yet been looked at within the context of the fragmentation model discussed 
above. It is therefore unclear how rapidly the HII region will grow, and whether it 
will immediately act to shut off accretion, or whether pockets of dense, cold gas can 
survive within the HII region for an extended period. 

• Do any low-mass Population III stars survive until the present day? 



One of the most exciting results of the Greif et al. ( 201 la] > model is that some of 



the protostars that are ejected from the centre of the star-forming minihalo have 
masses that are below 0.8 M , and hence lifetimes that are longer than the current 
age of the Universe. If these protostars avoided accreting any further gas, then they 
could have survived until the present-day, raising the possibility of directly detecting 
truly metal-free stars within the Milky Way. However, as we have already discussed 
above, the number of protostars with sub-solar masses that are ejected from the 
star-forming region is very sensitive to the way in which protostellar mergers are 
treated, and hence is highly uncertain at present. In addition, it is possible that any 
Pop. Ill protostars that have survived until the present day have also become too 
polluted with metals by ongoing accretion from the ISM for us to recognize them 



as metal-free stars, although the best current estimates (|Frebel, Johnson & Bromm 



2009 ; Johnson & Khochfar 201 l| i suggest that the effects of pollution will be very 



small. 



4 Summary 



In this review, we have focussed on three main topics: how the first star-forming 
minihalos come into existence and why they have the properties that they do; how 
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gas within a representative minihalo cools, collapses and forms a protostar; and how 
this protostar and the massive clump of gas surrounding it subsequently evolve. 

We have seen that on large scales, we now have a relatively good understanding 
of the physical processes involved in the formation of the first star-forming miniha- 
los. In order for gas to accumulate within a dark matter minihalo, it must be able to 
overcome the effects of both gas pressure and also the large-scale streaming motion 
of the gas relative to the dark matter. Since this streaming motion is typically su- 
personic, the latter effect generally dominates, and the result is that gas is prevented 
from accumulating in large quantities within minihalos with masses of less than 
around 1O 5 M . Within more massive minihalos, the gravitational force exerted by 
the dark matter is strong enough to overcome the effects of the gas pressure and the 
coherent streaming, and the gas begins to undergo gravitational collapse, reaching 
densities that are several hundred times higher than the cosmological background 
density. As the gas collapses, however, it is heated by compression and shocks. In 
order for the collapse to continue, the gas must be able to dissipate this energy, 
which it does through rotational and vibrational line emission from H2. 

In Section [L2| we saw that the amount of H2 formed within a given minihalo 
is a strong function of the temperature of the gas, with the final molecular fraction 
scaling roughly as xh 2 « r 3 / 2 with the temperature T. The H2 cooling rate is also a 
strong function of temperature. As a result, one finds that there is a critical minihalo 
virial temperature, r cr j t ~ 1000K marking the division between cooler halos that do 
not dissipate much energy within a Hubble time, and hence which do not form stars, 
and warmer halos that do manage to cool and form stars. As the virial temperature 
of a minihalo is a simple function of its mass and redshift, one can derive a critical 
minihalo mass that must be exceeded in order for the gas to cool effectively. This 
critical mass scales approximately as M crA ~ 1.6 x 10 6 ( 1 + z/10)~ 3//2 M , given 
standard values for the cosmological parameters. Combining this constraint with 
that arising from coherent streaming, one finds that at redshifts z > 40, the minimum 
mass of a star-forming minihalo is set by the need to overcome the effects of the 
streaming, and is roughly 10 M©, while at z < 40, H2 cooling is the limiting factor, 
and the minimum mass scale is somewhat larger. 

On smaller scales, we have also developed an increasingly good understanding 
of how the gas evolves as it cools, undergoes runaway gravitational collapse, and 
forms the first protostar. As outlined in Section[2] the gas first passes through a "loi- 
tering" phase, during which cold gas accumulates at the centre of the minihalo. The 
temperature and density of the gas at this point depend on the nature of the dom- 
inant coolant. When H2 dominates, we have T ~ 200 K and n ~ 10 4 cm~ 3 , while 
if HD dominates, then T ~ 100 K and n ~ 10 6 cm . The loitering phase ends and 
the collapse of the gas accelerates once the mass of cold gas that has accumulated 
exceeds the local value of the Bonnor-Ebert mass, which is around 1000 M Q in 
the H2-dominated case, but only 40 M© in the HD-dominated case. The next major 
event to occur is the onset of three-body H2 formation at n ~ 10 8 cm~ 3 which rapidly 
converts most of the atomic hydrogen into H2. The associated heat input leads to an 
increase in the gas temperature to T ~ 1000-2000 K, with the details depending to 
a significant extent on the rate coefficient chosen for reaction 41 which is poorly 
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constrained at low temperatures. At n <~ 10 10 cm~ 3 , the gas becomes optically thick 
in the main H2 cooling lines, but remains optically thin in the continuum. It can 
therefore continue to cool reasonably effectively at these densities, with the mean 
temperature only rising relatively slowly with increasing density. At n <~ 10 14 cm~ 3 , 
a new process, collision-induced emission from H2, begins to dominate the cooling. 
However, this does not lead to a significant drop in the gas temperature, as the gas 
quickly becomes optically thick in the continuum. At densities above n <~ 10 16 cm~ 3 , 
further radiative cooling of the gas is ineffective and the only remaining process ca- 
pable of slowing the temperature rise is collisional dissociation of the H2. While the 
H2 fraction in the gas remains significant, the temperature is prevented from rising 
much above 3000 K, but once most of the H2 has been destroyed, the temperature in 
the core rises steeply, and the internal thermal pressure eventually becomes strong 
enough to halt the collapse. State-of-the-art simulations have followed the gravita- 
tional collapse of the gas up to this point, which we can identify as the moment at 
which the first true Population III protostar forms. 

Nevertheless, several uncertainties remain in this picture of Pop. Ill star forma- 
tion. As already noted, the uncertainty in the three-body H2 formation rate limits 
the accuracy with which we can model the chemical and thermal evolution of the 
collapsing gas. In addition, current three-dimensional collapse models make use of 
simplified treatments of the effect of opacity on the H2 cooling rate, and the un- 
certainty that this introduces into the models has not yet been properly quantified. 
Further uncertainty comes from two additional issues which have only recently be- 
gun to be addressed: the role played by magnetic fields, and the influence of dark 
matter annihilation. Although the strength of any seed magnetic field existing prior 
to the assembly of the first star-forming minihalos is still poorly constrained, it now 
seems clear that the small-scale turbulent dynamo acting during the collapse of the 
gas will rapidly amplify even a very weak initial field up to a point at which it could 
potentially become dynamically significant. However, neither the final strength of 
the field nor its correlation length are well constrained at present, and without a bet- 
ter understanding of these values it is difficult to say to what extent the magnetic 
field will influence the details of the collapse. The role played by heating and ion- 
ization due to dark matter annihilation is even less well understood. Simple models 
suggest that it may be extremely important and may result in the formation of "dark 
stars" supported by the energy released by dark matter annihilation rather than by 
nuclear fusion, but the only hydrodynamical study performed to date suggests that 
the influence on the collapse is small, and that dark stars do not actually form. 

Finally, there remains the question of how the gas evolves after the formation 
of the first protostar. For much of the last decade, the leading model for this phase 
of the evolution of the gas has been what we have termed the "smooth accretion" 
model. In this model, it is assumed that the gas surrounding the newly formed pro- 
tostar does not fragment, but instead simply smoothly accretes onto the protostar, 
primarily via a protostellar accretion disk. Considerable work has been done within 
the framework of this model to understand the structure of the protostar during the 
accretion phase, and the effect of protostellar feedback on the surrounding gas. This 
work has shown that any feedback occurring prior to the protostar joining the main 
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sequence is unlikely to significantly reduce the accretion rate, and that the most 
plausible mechanism for terminating the accretion is photoionization of the accre- 
tion disk by ionizing radiation from the central star, implying that it must already 
have grown to some tens of solar masses. This model therefore predicts that Pop. Ill 
stars will generally be solitary, with only one or two forming in each minihalo, and 
massive, with masses M ^> 10M Q . 

Over the past couple of years, however, several new studies have appeared that 
have cast considerable doubt on the smooth accretion model. These studies have 
attempted to directly model the evolution of the gas as it begins to be accreted, and 
have shown that the accretion disk that builds up around the protostar is unstable 
to gravitational fragmentation even if the stabilizing effects of accretion luminosity 
feedback from the central protostar are taken into account. Once a few fragments 
have formed, the dynamical interactions between the individual fragments and be- 
tween the fragments and the gas can lead to further fragmentation, and to the ejection 
of low-mass fragments from the system. If we assume that all of the gravitationally 
bound fragments form protostars, then the result of this model is the assembly of a 
small, extremely dense cluster of Pop. Ill protostars with a wide range of masses. 



As discussed in Section 3.2.4 many aspects of the fragmentation scenario remain 
unclear and much work remains to be done before we can hope to have a good under- 
standing of the final protostellar mass function. Nevertheless, these results suggest 
that Population III star formation perhaps has far more in common with present-day 
star formation than has been previously recognised. 
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